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NOMENCLATURE 


Latin Symbols 

A reaction rate constant; also area, m 2 

C concentration, kg.mole/m 3 

C p specific heat, J/(kg.K) 

D diffusion coefficient, m 2 /s 

E total internal energy, J/kg; also activation energy, J/kg 

f mass fraction 

g Gibbs energy, J/(kg.K) 

h static enthalpy, J/kg 

h R base enthalpy, J/kg 


L; spectral radiative intensity, kW/(m 2 .sr.cm' 1 ) 

k thermal conductivity, J/(m.s.k); also line intensity to spacing ratio, 

cm* 1 , atm ' 1 

k b backward rate constant 

keq equilibrium constant 

kf forward rate constant 

L nozzle length, m 

L m mean beam length, m. 

m^ total number of narrow bands 

M molecular weight 

N temperature coefficient in reaction rate expression 

N s number of species 

N r number of reactions 

p gas pressure, atm 

flew conductive wall flux, kW/m 2 

—V.qr radiative source term, kW/m 3 


vii 


i 


Qrw 

net radiative wall flux, kW/m 2 

f- 

Q 

radiative energy per unit volume, kW/m 3 

r 

R 

gas constant, J/(kg.K); also random number 

i 

R u 

universal gas constant, J/(kg.K) 

r 

s, s',s" 

position variables, m 

h 

\ 

t 

time, s 

!! 

T 

absolute temperature, K 

1 

u 

velocity in x direction, m/s 

h. 

Um 

mean velocity, m/s 

If- 

u 

diffusion velocity in x direction, m/s 


V 

velocity in y direction, m/s 

1 

V 

diffusion velocity in y direction, m/s 


V 

diffusion velocity vector, m/s 

| 

w 

production rate of species, kg/(m 3 .s) 

r 

X 

x-coordinate, m 

i 

t 

X 

mole fraction 


y 

y-coordinate, m 

I , 

yb 

Greek symbols 

half height of cross sectional area of nozzle, m 

r 

t 

9 

line width to spacing ratio 


7 

stoichiometric coefficient; also half-width of an absorption line, cm -1 

t. 

£ 

equivalent line spacing, cm -1 


e 

polar angle 

[ 

p 

dynamic viscosity, kg/(m.s) 

r 

/Cp 

Planck mean absorption coefficient 

f 

t 


spectral absorption coefficient 


{,-0 

computational coordinates 

l 

P 

density, kg/m 3 


a 

normal stress, N/m 2 

r, 

k 

T 

shear stress, N/m 2 


TOJ 

spectral transmittance 

E 

<t> 

equivalence ratio 
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azimuthal angle 
wavenumber, cm 
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Chapter 1 
INTRODUCTION 


Extensive research is underway at the NASA Langley Research Center to develop 
hydrogen-fueled supersonic combustion ramjet (scramjet) propulsion systems for National 
Aero-Space Plane (NASP). A critical element in the design of scramjets is the detailed 
understanding of the complex flowfield present in the different regions of the engine over 
a wide range of operating conditions. Numerical modeling of the flow in various sections 
has proven to be a valuable tool for gaining insight into the nature of these flows [l^l]. 

In a hypersonic propulsion system, combustion takes place at supersonic speeds to re- 
duce deceleration energy losses. The products of hydrogen-air combustion are gases such 
as water vapor and hydroxyl radicals. These species are highly radiatively absorbing and 
emitting gases. Thus, numerical simulations must handle correctly radiation phenomena 
associated with supersonic flows. 

Over the past 30 years the analysis of radiative heat transfer has received increasing 
attention. This was first due to the advent of the space age, which made it necessary 
to develop tools to predict heat transfer rates in such high-temperature applications as 
rocket nozzles and space vehicle reentry, and in vacuum applications for spacecraft in 
outer space. Following a lull during the 1970s and early 1980s, interest in radiative heat 
transfer has recently increased again because of the need to predict and measure heat 
transfer rates in ever higher temperature applications in furnaces, and MHD generators, 
as well as the scramjet mentioned earlier [5-7]. 

Among the three modes of heat transfer, radiative heat transfer is quite different 
from conductive and convective heat transfer. Under normal conditions, conduction and 
convection are short-range phenomena. Thus we are able to perform an energy balance 
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in an infinitesimal volume. The principle of conservation of energy then leads to a 
partial differential equation. This equation may have up to four independent variables 
(three space coordinates and time). Thermal radiation, on the other hand, is generally a 
long-range phenomena [6-8]. Thus, conservation of energy cannot be applied over an 
infinitesimal volume, but must be applied over the entire volume under consideration. 
This leads to an integral equation involving up to seven independent variables (the 
frequency of radiation, three space coordinates, two coordinates describing the direction 
of travel of photons, and time). 

The analysis of thermal radiation is complicated further by the behavior of the 
radiative properties of materials. Properties relevant to conduction and convection are 
fairly easily measured and are generally well behaved. But radiative properties are 
usually difficult to measure and often display erratic behavior. For liquids and solids, the 
properties normally depend only on a very thin surface layer, which may vary strongly 
with surface preparation and often change from day to day. All radiative properties (in 
particular for gases) may vafy strongly with wavenumber, adding another dimension to 
the governing equation. Rarely, if ever, can this equation be assumed to be linear. 

Because of these difficulties inherent in the analysis of thermal radiation, accurate 
prediction of radiation in most realistic systems is currently still out of the question, 
although tremendous efforts have been made and significant progress has been achieved 
in the past decades. Prior to the 1970s, radiative transfer analyses were limited to one- 
dimensional formulations. Even for one-dimensional cases, nongray radiative heat trans- 
fer formulations were very complicated and their solutions required enormous amount of 
computational resources. Important works in nongray one-dimensional formulation have 
been reviewed in Refs. 5, 8-10. Since the 1970s, efforts have been directed toward 
formulating multi-dimensional equations for radiative transfer. Great achievements have 
been made for gray gaseous systems. However, studies on multi-dimensional nongray 
gaseous systems encounter tremendous difficulties and little progress has been made so 
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far. A survey of various methods for multi-dimensional radiative transfer analysis has 
been made by Howell [11, 12]. Discussions were made regarding the feasibility of in- 
corporating spectral integration in the techniques using narrow band [13] and wide band 
models [14, 15]. Another review [16] has provided details of several methods that could 
possibly be applied to multi-dimensional radiative transfer in molecular participating me- 
dia. Different review articles have indicated unanimously that one of the most promising 
methods to investigate nongray participating media in multi-dimensional systems is the 
Monte Carlo method (MCM). 

The MCM is a statistical sampling technique which can simulate exactly all impor- 
tant physical processes. In this method, the numerical treatment of the mathematical 
formulation is easy and the usual difficulties encountered in complex geometries can be 
circumvented easily. It is because of these advantages that the MCM has been applied 
to solve many radiative transfer problems. The earliest application of this method for 
radiative transfer problems was made by Howell and Perlmutter [17]. Radiative problems 
of increasing complexity which have been investigated by this method have appeared in 
the literature [18-22]. Studies on reducing the computational time by using this method 
are also available [23, 24]. The gray gas assumption, however, is made in most of these 
analyses. 

Like other numerical methods, the MCM also has some disadvantages. One of them 
is the large appetite for computer time, and another is the statistical fluctuation of the 
results. With the rapid development of computers, these two disadvantages are becoming 
of less concern and interest in the MCM is increasing. One of the recent applications of 
the MCM has been in the investigation of radiative interactions in nongray participating 
media using a narrow band model. For example, Taniguchi et al. [25] applied a simplified 
from of the Elsasser narrow band model to investigate the problem of radiative equilibrium 
in a parallel plate system. Farmer and Howell [26] obtained a Monte Carlo solution of 
radiative heat transfer in a three-dimensional enclosure with an anisotropically scattering, 
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spectrally dependent, inhomogeneous medium. Modest [27] discussed the effects of 
narrow band averaging on surface and media emissions. It was pointed out that the 
narrow band model may be applied successfully to the MCM after verification in an 
isothermal and homogeneous medium. However, all these studies have failed to reflect 
some fundamental mechanisms of the MCM in conjunction with a narrow band model, 
and the application of the MCM to nongray radiation problems is still uncertain. 

The first objective of this study is to employ a general and accurate narrow band model 
to investigate radiative heat transfer using the MCM. The same nongray model has been 
applied to investigate radiation contributions using the discrete direction method [28] and 
the S-N discrete ordinates method [29]. The present investigation includes derivation 
of the Monte Carlo statistical relationships, discussion of the fundamental features that 
are different from other methods and demonstration of the capability of the MCM for 
nongray analyses. A one-dimensional problem is considered first, and the validation 
of the Monte Carlo analysis is conducted by comparing the Monte Carlo results with 
available solutions for the cases with and without other modes of heat transfer. Next, 
the Monte Carlo formulations suitable for multi-dimensional problems are developed 
and validated. From our knowledge, the present study is the first to provide accurate 
and general radiative transfer formulations which are applicable for any nongray and 
multi-dimensional system. 

A literature survey indicates that a great deal of effort has been made towards 
an accurate formulation of radiative transfer equations. Most applications of these 
formulations have been restricted to non-reacting homogeneous systems. Only a limited 
number of studies are available to investigate the interaction of radiation heat transfer 
in chemically reacting, viscous, compressible flows such as those in scramjet propulsion 
systems. Mani and Tiwari [30] were the first to take into account the effect of radiation 
on chemically reacting supersonic flows. This work has been extended to include some 
relatively more advanced chemistry models by Tiwari et al. [31]. In both of these studies. 
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a tangent slab approximation for radiative transfer was employed. This approximation 
treats the gas layer as a one-dimensional slab in evaluation of radiative flux. Obviously, 
it is impossible to obtain reliable quantitative predictions of radiative heat transfer from 
this treatment. Therefore, the second objective of this study is to apply the Monte 
Carlo formulations developed during the course of present research efforts to investigate 
the radiative interaction in multi-dimensional chemically reacting flows. The specific 
problem considered is the supersonic flow of premixed hydrogen and air in an expanding 
two-dimensional nozzle. Two-dimensional radiative heat transfer in this problem is 
simulated using the MCM; the results of radiative flux are then incorporated in the two- 
dimensional Navier-Stokes equations. This procedure provides a more accurate prediction 
of radiative effects on flowfield and wall heat transfer than those available in previous 
studies. The physics of radiative interactions in chemically reacting compressible flows 
can be understood more clearly from this study. 

Two different objectives divide the present study into two parts. The first part is 
to develop and validate the Monte Carlo formulations with a narrow band model. This 
work is included in Chaps. 2-5. Information on radiation absorption models is given 
in Chap. 2. Development and validation of the Monte Carlo formulations for one- 
dimensional problem is provided in Chap. 3. Further validation for the one-dimensional 
formulations is conducted in Chap. 4 by considering a simple problem of radiative 
interactions. Development and validation of Monte Carlo formulations for radiative 
transfer in multi-dimensional systems is presented in Chap. 5. The second part of the 
study is an investigation of radiative interactions in chemically reacting flows. This work 
is included in Chap. 6. Finally, the conclusions reached from this study are summarized 
in Chap. 7. 
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Chapter 2 

RADIATION ABSORPTION MODELS 

The study of radiative transmission in nonisothermal and inhomogeneous gaseous 
systems requires a detailed knowledge of the absorption, emission and scattering charac- 
teristics of the specific species under investigation. In absorbing and emitting media, an 
accurate model for the spectral absorption coefficient is of vital importance in the correct 
formulation of the radiative flux equations. A systematic representation of the absorption 
by a gas, in the infrared part of spectrum, requires the identification of the major infrared 
bands and the evaluation of the line parameters (line intensity, line half-width, and spac- 
ing between the lines) of these bands. The line parameters depend upon the temperature, 
pressure and concentration of the absorbing molecules and, in general, these quantities 
vary continuously along a nonisothermal and inhomogeneous path in the medium. In 
recent years, considerable efforts have been expended in obtaining the line parameters 
and absorption coefficients of important atomic and molecular species [32-34]. 

For an accurate evaluation of the transmittance ( or absorptance) of a molecular 
band, a convenient line model is used to represent the variation of the spectral absorption 
coefficient. The line models usually employed are Lorentz, Doppler, and Voight line 
profiles. A complete formulation ( and comparison) of the transmittance and absorption 
by these line profiles has been given [9, 10, 35-37]. In a particular band consisting of 
many lines, the absorption coefficient varies rapidly with frequency. Thus, it becomes 
a very difficult and time-consuming task to evaluate the total band absorption over the 
actual band contour, by employing an appropriate line profile model. Consequently, 
several approximate band models (narrow as well as wide) have been proposed which 
represent absorption from an actual band with reasonable accuracy [9, 10, 35-46]. Several 
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continuous correlations for total band absorption are available in the literature [9, 10, 
35-37, 43-46]. These have been employed in many radiative transfer analyses with 
varying degree of success [9, 10, 35-37, 47]. A brief discussion is presented here on 
narrow band models, wide band models, and band absorptance correlations. 

The absorption within a narrow spectral interval of a vibration-rotation band can be 
represented quite accurately by the so-called “narrow band models.” The most commonly 
employed narrow band models are the Elsasser, statistical, random-Elsasser and quasi- 
random narrow band models. Various narrow band models have been tested with the 
results of line-by-line calculations in the literature [37, 48, 49]. Accurate results for 
temperature and heat flux distributions were obtained with the statistical narrow band 
model, which assumes the absorption lines to be placed randomly and the intensities to 
obey an exponential-tailed-inverse distribution. The transmittance predicted by this model 
in a homogeneous and isothermal column of length 1 due to gas species j, averaged over 
[uj — (Aw/2), w+(Acj/ 2)], is expressed as [50] 


rl = exp 




2x Xjplk\ 



( 2 . 1 ) 


where Xj represents the mole fraction of the absorbing species j and p is total pressure; 
k and (3 — 2 ^ 7 /S are the band model parameters which account for the spectral structure 
of the gas. The overbar symbol indicates that the quantity is averaged over a finite 
wavenumber interval Aw. The narrow band width considered is usually 25 cm' 1 . 
Parameters k and 1/8 generated from a line-by-line calculation have been published 
for H 2 O, CO 2 , CO, OH, NO, and other species [33, 48, 51]. The mean half-widths 7 for 
H 2 O and CO 2 are obtained by Soufiani et al. as [48] 


1H 2 o = 0.066— il.OXHiOTf? + [1 -2(Xh 2 o + Xn 2 ) 
Ps ( 7 

fTrX) 


+0.8Xo 2 + 1 ,6Xco 2 


( 2 . 2 ) 
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and 

lco 2 = l®-^Xco 2 + 0.058(X/y 2 + Xo 2 ) + 0.15 Xh 2 o] (2.3) 

where p s and T s designate standard pressure and temperature(l atm, 296 K). Alternative 
formulations for evaluating the mean half- width 7 are also available in Ref. 33. 

The absorption within the spectral range of the entire vibration-rotation band can be 
represented by the so-called “wide band models.” The total band absorption of the wide 
band models is given by 

00 

A— J [1 — exp{— tt u l)]d(u> — u?o) (2.4) 

—00 

where the limits of integration are over the entire band pass, is the spectral absorption 
coefficient, and loq is the wave number at the center of the wide band. 

Four commonly used wide band models are the box, modified box, exponential and 
axial wide band models. The exponential wide band model, first developed by Edwards 
and Menard [41], is by far the most successful of the wide band models. In this model, 
the line intensity is assumed to be an exponential decaying function of the wave number 
[9, 10, 35-37, 44-46], such that 

§j_ _ _S_ e ~b 0 \u;-u;o\/Ao (2 5 ) 

d Aq 

where S is the band intensity, Ao the band width parameter and bo=2 for a symmetrical 
band or bo=l for bands with upper and lower wave number heads at uq . 

The radiative flux term usually involves multiple integrals even for simple geometries. 
As a result, numerical calculation of radiative flux for energy transfer becomes very time 
consuming. Therefore it is desirable to replace the relation for the total band absorptance, 
given by Eq. (2.4), with a continuous correlation [5, 10, 52]. Numerous correlations are 
available in the literature for wide band absorptance. The first correlation to satisfy the 
linear, square-root, and logarithmic limits of the wide band absorptance was proposed 
by Edwards and Menard [41]. The most widely used correlation is the Tien and Lowder 
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continuous correlation because of its simplicity and relative accuracy. This correlation 
is expressed as [52] 


A = Aq In 




(u 4- 2) 

(u + 2 /(<)) 



( 2 . 6 ) 


where 

f(t) = 2.94[1 - exp(-2M)], t = y (2.7) 

Here u=Spl/Ao is nondimensional path length and /?* is line structure parameter. Wide 
band model correlation parameters for various gases are available in the literature [6, 
7, 10, 35]. 

Among the approximate band models discussed above, the wide band models and 
band absorptance correlations are simpler than the narrow band models, and they have 
been used extensively in the study of nongray radiative heat transfer for the past three 
decades. However, the spectral discretization used in the wide band models and band 
absorptance correlations is too wide and it does not take into account the low resolution 
correlations between intensities and transmissivities. This leads to significant temperature 
and heat flux discrepancies [49]. Also, the case of partially reflecting walls cannot be 
modelled correctly with these models [10]. Recently, the narrow band models have begun 
to receive attention due to the rapid development of computers and strong requirements 
for accurate analyses of radiation [25-29]. Some narrow band models compare favorably 
to the line-by-line calculations; However, they are much simpler than the line-by- 
line models. In addition, use of the narrow band models can avoid some notorious 
disadvantages occurring with the wide band models and band absorptance correlations. 
In this study, the narrow band model expressed in Eq. (2.1) is employed to investigate 
nongray radiative heat transfer. 

For a nonisothermal and inhomogeneous column, the Curtis-Godson approximation 
[53] leads to accurate results if pressure gradients are not too large. Basically, this 
approach consists of transformation of such a column into an equivalent isothermal and 
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homogeneous one. For the narrow band model expressed in Eq. (2.1), effective band 
model parameters k e and (3 e are introduced by averaging k and (3 over the optical path 
U of the column as 

i 

U(l) = J p(y)X j (y)dy (2.8) 

0 

l 

k e = jjfic J p{y)Xj(y)k(y)dy (2.9) 

0 

l 

fie = J p(y)X j (y)k(y)fi(y)dy (2.10) 

The transmittance of this equivalent column is then calculated from Eq. (2.1). 

A distinguishing characteristic for the band models discussed above is the dependence 
of the wavenumber. If it is assumed that the absorption coefficient is independent of the 
wavenumber, the radiation absorption is then represented by the so-called “gray model”. 
The gray model is rarely a physically realistic approximation, but it serves as an initial 
step for studying the effect of radiative heat transfer. For a nonuniform temperature field, 
the gray model used for optically thin radiation is the modified Planck mean absorption 
coefficient which, for black bounding surfaces, is defined as [8, 35] 

Km(T,T w ) = k p (T w )(T w /T) (2.11) 

where /c p (T) represents the Planck mean absorption coefficient. For a multiband system 
of a homogeneous gas, /c p (T) is expressed as 

n 

k p (T) = [e t (a, i ,T)5,(T)]/(ar 4 ) (2.12) 

J— 1 

where n represents the number of vibration-rotation bands, eb(tJi, T) is the Planck function 
evaluated at the ith band center, Sj(T) is the integrated band intensity of the ith band, and 
a is the Stefan-Boltzmann constant. Equation (2.12) is modified to apply to a mixture 
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of different gases as 

k p (T) = j£ [e»(c*,T)Si(r)] j /(aT*) ( 2 . 13 ) 

where j denotes the number of species in the mixture and pj is the partial pressure of 
jth species. The band model parameters for various gases are available in the literature 
[ 6 , 7 , 10 , 35 ]. 
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Chapter 3 

MONTE CARLO SIMULATION USING A NARROW BAND MODEL 

In this chapter, the radiative heat transfer for a one-dimensional problem is investi- 
gated using a narrow band model and Monte Carlo simulation. The physical model is 
established in Sec. 3.1. Monte Carlo formulations are developed in Sec. 3.2. Discussion 
of special features of the method for nongray analyses is made in Sec. 3.3. Estimation 
of statistical error is presented in Sec. 3.4 and validation of Monte Carlo formulation 
is conducted in Sec. 3.5. 


3.1 Physical Model 

To investigate radiative heat transfer using the MCM and a narrow band model, 
a simple problem is considered at first Figure 3.1 shows an absorbing and emitting 
molecular gas between two infinite parallel plates with slab thickness of L. Temperature, 
concentration and pressure in the medium are assumed known. The walls are assumed to 
be diffuse but not necessarily gray. The wall temperature is also assumed known. Usually, 
the radiative transfer quantities of interest are the radiative source term — V.g r inside the 
medium and the net radiative wall flux q TW • In order to calculate these quantities, the 
medium considered is divided into (M— 2) volume elements. The grid 1 and M are 
numbered on the lower and upper walls, respectively. Temperature, concentration and 
pressure are assumed to be constant in each volume element. The typical method for 
handling radiative exchange between surface and/or volume elements is to evaluate the 
multiple integral, which describes the exchange, by some type of numerical integration 
technique. This usually is a good approach for simple problems, but an alternate method 
is used here. Radiative transfer in the computational domain is simulated using the MCM. 
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For an arbitrarily chosen volume element with a volume, 8V> and an arbitrarily chosen 
surface element with an area 8 A, the relations for — V.g r and q rw are expressed as 


-V.tfr 


Qv-sv + Qa-bv - Qsv 
8V 


(3.1) 


Qv-SA + Qa-SA - QsA ^ 

qrw = — (3.2) 

Here, Qv-sv and Qvsa are the total radiant energy from the entire gas that is absorbed 
by the volume element 8V and surface element 8 A, respectively; Qa-sv and Qa-sa 
are the total radiant energy from the bounding walls that is absorbed by 8V and 8 A, 
respectively; Qsv and Qsa are the radiant energy emitted by 8V and 8 A, respectively. 

To evaluate the terms Qv-sv > Qa-sv* Qsv and Qv-sa in Eqs. (3.1) and (3.2), 
the MCM uses a large number of bundles of energy (statistical samples) to simulate 
the actual physical processes of radiant emission and absorption of the energy occurring 
in the medium. These energy bundles are similar to photons in their behavior. The 
histories of these energy bundles are traced from their point of emission to their point of 
absorption. What happens to each of these bundles depends on the emissive, scattering 
and absorptive behavior within the medium which is described by a set of statistical 
relationships. The total number of energy bundles absorbed by each element multiplied 
by the energy per bundle gives the interchange of radiation among the volume and /or 
surface elements. The values of — V.<? r and q rw can then be obtained from Eqs. (3.1) 
and (3.2), respectively. 

The use of a narrow band model in the MCM presents new features in the analysis of 
radiative heat transfer. The statistical relationships currently in use need to be modified. 
The following Monte Carlo analysis is based on an arbitrarily chosen finite volume 
element. The statistical relationships for an energy bundle emitted from a surface element 
can be derived by following the same procedure. 
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3.2 Monte Carlo Formulations 


Let us consider the Planck spectral blackbody intensity Ibo> that enters the ith volume 
element at the point s on the lower side and intersects the upper side at the point s' as 
shown in Fig. 3.1. A spherical coordinate system is established and centered at the point 
s. From Ref. 6, the amount of energy emitted for a wavenumber range du and along a 
pencil of column s— >s' with a solid angle increment dO is expressed as 

dQi = hu> [l ~ t w (s — * s 1 )] cos QdSldw (3.3) 


where rm(s— ►s') is the spectral transmittance over the path s— >s', 0 is the polar angle 
between the y axis and the direction of the column s— >s', and dft=sin0d0dV> where ip is 
the azimuthal angle. The total emitted energy per unit volume is obtained by integrating 
Eq. (3.3) over all wavenumbers, and polar and azimuthal angles as [6, 54, 55] 


oo r 2tt 


«■-/// hu [l — (s — > s ; )] cos 6 sin OdipdOduj 
ooo 

OO IT 

= 27T //'-[■ — r w (s — ► s')] cos 0 sin OdOdu 


o o 

oo 1 


= 27 T J J ~ T u (Ayi/fi)]fJidfJt<kj; fi = cos 6 (3.4) 


o -1 


where Ay* is the thickness of ith volume element. It should be noted that the sign of 
Ayi is different when fi varies from positive to negative. 

The simulation of an energy bundle includes the determination of wavenumber and 
direction of emission of this energy bundle in the finite volume element. The statistical 
relationships for determining these parameters are readily obtained from Eq. (3.4) as 


Rw = 


U) 1 

2 t rff I bu) [ 1 - T u (Ayi/ fj,)]fidfj,duj 
o -l 

Qi 


(3.5) 


1 oo 


R. 


2 vf J Iiul 1 - T w (Ayi/fi)]fidudfi 
» o 

Qi 


(3.6) 
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where Ro; and are random numbers which are distributed uniformly between zero 
and one. In Eqs. (3.4)-(3.6), r w is a real spectral transmittance. Before solving these 
equations to obtain u and fj, from a set of given values of Ru; and R^, the narrow band 
model should be applied to approximate the real spectral transmittance. 

For the narrow band model, the absorption bands of the gas are divided into spectral 
ranges Au> wide; each is centered at and characterized by the superscript k; the band 
parameters obtained are the averaged quantities over a narrow band. So, the spectral 
quantities in Eqs. (3.4)-(3.6) should be transformed into the averaged quantities over 
a narrow band for practical applications. Taking the spectral average over all narrow 
bands, Eqs. (3.4)-(3.6) are expressed as 


rtiu, i \ 

Qi — 27r ^ < / I b(J k[l -T^(Ayi/fi))fidfj, 
k=i t_i 


AxJ 


(3.7) 


X) ( f hu>*{^ ~ ^{Ayi/ n)\iidn Iaw* 

Ru, = — - l— tt , (w- 1 <u<u n ) (3.8) 

V* 


TTt(jj 


Ru — 


2tt E { MuM 1 - r w k{Ayi/ fi)\fxdfx \Au k 
k=l 

Qi 


(3.9) 


where m^ is the total number of narrow bands. The following narrow band approximation 
has been used in obtaining Eqs. (3.7)-(3.9) 

^bui k ^~w k j dtO 

Au k 


hu k 


(a w * / 

\ Au> k 


r L ,dio 


= hu kT w k (3*10) 

This is because Ibu> is essentially constant over a narrow band and may be taken out of 


the spectral integral. Otherwise, the average product hu kT u k is not equal to the product 


of I buJ k and T^fc. 
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Equations (3.8) and (3.9) are solved for lj and fi each time a set of values of Ru> and 
are chosen. The computing time becomes too large for practical calculations since the 
integrands in these equations are very complex functions of integration variables and the 
number of energy bundles usually is very large. To circumvent this problem, interpolation 
and approximation methods are employed. For example, to obtain the value of lj for a 
given value of R^>, we first choose different values of uj and obtain the corresponding 
values of Ru> from Eq. (3.8). Then, a smooth curve is constructed to match these data 
points, and lj values are obtained easily from this curve for selected values of R^. The 
procedures for determining // are similar to those for lj. 

Following the determination of wavenumber and direction of an energy bundle, it 
is essential to find the location of absorption of the energy bundle in the participating 
medium. Let us still consider the emitted radiant energy along a pencil of column 
s— >s' (Fig. 3.1). After this amount of energy is transmitted over a column s'— ►s", the 
remaining radiant energy is given by 

dQ[ = hu>[ 1 — t u (s — ► <s / )]t w ( < s / — ► s") cos Odftdcj (3.11) 

where t w (s'— ►s") is the spectral transmittance over the path s'— >s". Taking a narrow 
band average over Eqs. (3.3) and (3.11) and dividing the latter one with the first one, the 
statistical relationship for determining the location of absorption can be expressed as 

K - i 1 ~ T »( 8 3 ") 

1 1 - t^( s -> s') 

_ gjT(V -» S'') ~ Tu(s -» S^T^S* S") 

1 ( * ^ 

where R\ is a random number. The averaged product r w (s — ► s 1 )^ (s' — > s") is not equal 
to the product of —*■ s 1 ) and s 1 — > s") because the t u (s —+ s 1 ) and r a ,(5 , — > s n ) 

have a strong wavenumber dependence due to the high resolution structure in a very 
small range of an absorption band (hundreds of major absorption lines in a 25 cm -1 
spectral interval), and must be treated in a spectrally correlated way. Equation (3.12) 
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can be simplified as 


/? _ T( *( s ' s ") ~ T “( s s ") 
l ~ 

If the spectral correlation between 7^(3 — > s') and t^s' — > 3") 
then Eq. (3.12) becomes 


(3.13) 

is not taken into account. 


Rl = r a ,( 3 / — > s") 


(3.14) 


Equation (3.14) is the statistical relationship usually employed for determining the 
location of absorption in the Monte Carlo simulation and is quite different from Eq. 
(3.13). For an isothermal and homogeneous medium, the travelling distance of an energy 
bundle can be obtained directly by solving Eq. (3.13) for a given random number. But this 
procedure turns out to be somewhat complicated for a nonisothermal and inhomogeneous 
medium. It becomes necessary to try each volume element starting from the adjacent 
element of the location where an energy bundle emits until a finite volume element is 
found in which Eq. (3.13) can be satisfied. 


33 Special Features of MCM for Nongray Analysis 

The MCM is quite different from other numerical techniques for the analysis of 
radiative heat transfer. Its characteristics have been discussed in detail by Siegel and 
Howell [6]. Use of a nongray model in the radiative transfer analysis requires significant 
changes. Two special features of incorporating the nongray model in the MCM will be 
discussed. 

Most of the existing analyses in radiative heat transfer start with the transfer equation 
of the type given by Siegel and Howell [6]. In order to apply a narrow band model, 
this equation has to be spectrally averaged over a narrow band. This averaging treatment 
results in two types of spectral correlations [56]. One is the spectral correlation between 
the intensity and the transmittance within the medium. Another is the spectral correlation 
between the reflected component of the wall radiosity and the transmittance. In order 
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to investigate the first type of spectral correlation, all the intermediate transmittances in 
each finite volume element of medium along the path the radiative energy travels must 
be calculated and stored to make a correlated calculation. In order to investigate the 
second type of spectral correlation, a series expansion of the wall radiosity is required 
[57, 58]. Essentially, this series expansion is utilized along with a technique for closure 
of the series. 

The simulation of radiative heat transfer in the MCM is not based directly on the 
radiative transfer equation. This results in the MCM having features different from the 
other methods for nongray analysis. When radiative energy is transmitted in a medium, 
the spectral correlation does occur in the MCM, but it occurs between the transmittances 
of two different segments of the same path which is different from other methods. This 
is the first noteworthy feature of the MCM for nongray analysis. 

The MCM procedures are based on the direct simulation of the path of an energy 
bundle. For the case with reflecting walls, the mechanism of the reflections simulation in 
the MCM is the same as a series expansion of the wall radiosity. However, this simulation 
process becomes much simpler because of its probabilistic treatment. Also, there are no 
spectrally correlated quantities involved. This is the second distinctive feature of the 
MCM for nongray analysis. Exact treatment of the reflections in the MCM in nongray 
gases is the same as that in gray gases and may be found in the literature [6, 54, 55]. 

The second feature of the MCM allows one to obtain results for a reflecting wall with 
very little increase in the computation time compared to that for a nonreflecting wall. But 
in other methods, the consideration of the history of a finite number of reflections and 
approximating the remaining reflections by a closure method in the radiative transfer 
equation complicates the mathematical formulation and increases the computer time 
considerably. As the geometry considered becomes more complicated, exact simulation 
of the radiative heat transfer in cases with reflecting walls will be very difficult for most 
existing methods, while it is not a big problem for the MCM. So, it seems that the MCM 
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is able to retain the feature of simplicity in dealing with complicated problems while a 
narrow band model is employed. 

3.4 Estimation of Statistical Error 

In the Monte Carlo simulation, the computational error consists of the statistical error 
and the computer truncation error. The statistical error is the major error source and the 
truncation error is usually neglected. From probability theory [59-63], the convergent 
speed of the Monte Carlo solution is proportional to the 1 /y/N for a statistical process 
with a sample size of N. Such a speed is very slow among all kinds of numerical 
computation methods. In practical applications, sample size cannot be an infinitely large 
number due to limitations on computer resources. Therefore, the Monte Carlo calculation 
must be supplemented with an estimate of the statistical error. 

To analyze the statistical error, the radiation simulation between two elements is 
considered first. For the sake of simplicity, the element from which the radiant energy is 
emitted is represented by SVi and the element from which the radiant energy is absorbed 
is represented by SVj ; it does not matter whether the element considered is a volume 
element or a surface element. 

In the computational domain, the travel state of an energy bundle emitted from 8V{ 
can be described by the spatial position r and moving direction Thus, the travel state 
is expressed as 

S = (r, n) (3.15) 

After travelling the i-th step, the state becomes 

5 = (rj, fii) (3.16) 

An energy bundle travels in the medium surrounded by the surfaces and is absorbed by 
an element 8Vj at the r-th step. Such a random travelling process X can be described 
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Thus, the statistical simulation error in the radiative term is expressed as 


€ = \¥V-6Vj - Qv-SVj\ < Xa^V-SVj (3.29) 

where x a is the confidence coefficient and a is the confidence level. Table 3.1 is 
the standard normal distribution table which provides the relation between the confidence 
coefficient xa and confidence probability 1-a. For example, when l-a=0.95, xa is taken 
to be 1.96. In practical applications, the relative statistical error is usually employed and 
it is expressed as 


£ £ Xa^V-BVj - 

Qv-svj K Qv-eVj ° 

where <$o is the maximum relative statistical error. 


(3.30) 


3.5 Results and Discussion 

In order to validate the Monte Carlo simulation, along with a narrow band model, 
results for a radiative source inside the medium and the net radiative wall heat flux have 
been obtained for different temperature and concentration profiles with nonreflecting and 
reflecting walls. Appendix A provides the computer code for the Monte Carlo simulation. 
In the present study, the reflectivities of two parallel diffuse walls are assumed to be 
identical and are denoted by the symbol p. Three different temperature profiles were 
used here: uniform, boundary layer type and parabolic profiles (Fig. 3.2). They were 
obtained from Kim et al. [29] and Menart et al. [56]. For the uniform temperature profile, 
the gas temperature was chosen to be 1000 K, while the walls were held at 0 K. Also 
shown in the figure is a parabolic H 2 O concentration profile for a mixture of H 2 O and N 2 
at 1 atm, and it was also taken from the above cited references. A uniform composition 
of pure H 2 O vapor at 1 atm is another H 2 O concentration profile that was used. Several 
cases with the selected temperature and H 2 O concentration profiles have been considered 
previously using the S-N discrete ordinates method by including all important bands 
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Table 3.1 Relation between confidence coefficient xcv and confidence probability 1 -a 


Xa 

1 - 0 ' 

Xcx 

1 -a 

Xa 

1 -a 

Xa 

1 -a 

0.0 

0.00000 

1.0 
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[29, 56]. The Monte Carlo solutions have been compared with published solutions for 
identical conditions. 

In the Monte Carlo simulation, the entire slab of the physical problem is divided 
into 20 sublayers for all calculations. Further subdivision of the computational domain 
was found to yield little change in the results. The computations were performed on a 
Sun Sparc workstation. The number of total energy bundles for each case was chosen 
to be 50,000. This choice represents a compromise between accuracy and economy of 
computation time. When the relative statistical errors of the results were chosen to be less 
than +3%, the probability of the results lying within these limits was greater than 95%. 
The computing times for the correlated and noncorrelated formulations were essentially 
the same. For an isothermal and homogeneous medium, the required CPU time was 
about 1-2 minutes for each case. For a nonisothermal and inhomogeneous medium, the 
CPU time increased to 5-7 minutes, and was nearly 10 minutes for the case with strongly 
reflecting walls (/9=0.9) with large optical length (L=0.5 m). 

The situation with nonreflecting walls is considered first. Figures 3.3-3.6 show the 
comparisons between the Monte Carlo solutions and S-N discrete ordinates solutions. 
Four different S-N discrete ordinates solutions are available in the literature [29] which 
employ different band models. For our comparison, we selected the solution — S-20 
nongray narrow band solution because it employs the same narrow band model as used 
in this study. 

Figures 3.3 and 3.4 show the radiative source results obtained for the uniform 
temperature and uniform pure H 2 O vapor distributions with slab thicknesses of 0. 1 m and 
1.0 m, respectively. The Monte Carlo results essentially match the S-N discrete ordinates 
results. Figure 3.5 presents the results for the boundary layer type temperature profile 
and for the same concentration distribution as in Figs. 3.3 and 3.4. The Monte Carlo 
results predict the same changes in gas behavior (from a net emitter near the hot wall to 
a net absorber away from the hot wall) as the S-N discrete ordinates results. The results 
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Fig. 3.6 Comparison of radiative source term for the parabolic H 2 O concentration profile. 
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for the parabolic H 2 O concentration distribution (with a uniform temperature profile) are 
shown in Fig. 3.6. The Monte Carlo method also predicts the interesting W type shape 
distribution of — V.<? r as in the S-N discrete ordinates method. Here the Monte Carlo 
solutions appear to be a little higher than the S-N discrete ordinates solutions, especially 
in the central region. 

The results for the net radiative wall heat flux obtained for the cases presented in 
Figs. 3. 3-3.6 are given in Table 3.2. The differences in results between the different 
solutions for the three cases are less than 3.5%. This shows agreement similar to that 
for the radiative source results. 

The situation with reflecting walls is considered next. Figures 3.7-3.11 show the 
comparisons between the Monte Carlo solutions and the S-N discrete ordinates solutions 
for different wall reflectivities and slab thicknesses. For these results, the parabolic 
type temperature profile and the uniform composition of pure H 2 O vapor at 1 atm were 
assumed. The S-N discrete ordinates solutions were based on the second-degree closure 
results [56]. The second-degree closure means that the history of two reflections is 
considered in the radiative flux equation and the remaining reflections are approximated 
by a closure method. Based on the study by Kim et al. [64], the second-degree 
discrete ordinates solutions for typical cases required about 160 minutes on a Cray-2 
supercomputer. This is significantly higher than the CPU time required for the MCM, 
which is not more than 10 minutes on a Sun Sparc workstation. 

Figures 3.7-3.9 present the results of —V.q r for the wall reflectivities of p= 0.1, 
0.5 and 0.9 respectively, with a slab thickness of L=0.5m. Excellent agreement between 
different solutions is seen in the figures. In the central region, the values of — V.# r are 
approximately constant. The Monte Carlo results appear to oscillate in that region. The 
reason is that the total number of energy bundles is a finite number and the Monte Carlo 
results are of a statistical nature. The oscillation decreases and the results of —V.q r 
become smoother as the total number of energy bundles is increased. These oscillations 
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Table 3.2 Comparison of net radiative wall heat fluxes with nonreflecting walls (kW/m 2 ) 



Monte Carlo 

S-N Discrete Ordinates 

Uniform T; L=0. 1 m 

-14.2 

-14.3 

Unifrom T; L=1.0 m 

-27.6 

-28.2 

Boundary layer T 

280.4 

277.4 

Uniform T with 

-24.5 

-25.4 

concentration profile 
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Fig. 3.7 Comparison of radiative source term in pure H 2 O for p= 0.1, L=0.5 m. 




Fig. 3.8 Comparison of radiative source term in pure H 2 O 
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Fig. 3.9 Comparison of radiative source term in pure H 2 O for p- 0.9, L=0.5 m. 
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Fig. 3.10 Comparison of radiative source term in pure H 2 O for p= 0.9, L=0.1 m. 
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are also exhibited in other figures. Figures 3.10 and 3.11 show the results for the strongly 
reflecting walls where p= 0.9, with slab thicknesses of L=0.1 m and L=1.0m, respectively. 
Again, the Monte Carlo solutions are in good agreement with the S-N discrete ordinates 
solutions. 

Table 3.3 shows the net radiative wall heat fluxes for the cases presented in Figs. 
3.7-3.11. The Monte Carlo results are slightly lower than the S-N discrete ordinates 
results. But the differences are within 6%. There are physical justifications for such 
discrepancies. In the S-N discrete ordinates method, the history of two reflections is 
taken into account and the remaining reflections are approximated as travelling in a 
medium without any attenuation. This approximation overpredicts the radiative energy 
absorbed on the walls. In the MCM, the history of the reflections is simulated in an exact 
manner. The Monte Carlo solutions are also subject to small statistical errors. 

The spectrally correlated results are compared with the noncorrelated results in 
Figs. 3.12 and 3.13. A spectral correlation has been considered in all the results 
presented in previous figures. In a spectrally noncorrelated formulation, the correlation 
between spectrally dependent quantities is neglected. By using Eq. (3.14), the Monte 
Carlo noncorrelated results can be obtained. The temperature and H 2 O concentration 
distributions considered here are the same as those in Figs. 3.7-3.11. The wall 
reflectivities are ^=0.0 for Fig. 3.12 and p=0.5 for Fig. 3.13, and slab thickness 
L is 0.1 m for both cases. The figures show clearly that the noncorrelated results 
overestimate the gas emission in the central region, and differ by about 30-35% from 
the correlated results. The reason for these discrepancies is in the derivation of the 
statistical relationship for determining the location of absorption of an energy bundle. 
The term rj[s -*■ s')r w («s' — ► s") in Eq. (3.12) can be treated in two different ways, that 
is, T w (s -> -> s") = nJ(s -+ s") and t u (s -» s') • r^s' -> s"), respectively. 

The first choice results in the correlated formulation given by Eq. (3.13) and the second 
choice results in the noncorrelated formulation given by Eq. (3.14). Since the value of 
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Table 3.3 Comparison of net radiative wall heat fluxes with reflecting walls (kW/m 2 ) 



L(m) 

Monte Carlo 

S-N Discrete 




Ordinates 

p=0A 

0.5 

14.42 

15.12 

II 

p 

0.5 

9.47 

9.66 


0.1 

2.22 

2.34 

p= 0.9 

0.5 

2.55 

2.70 


1.0 

2.58 

2.67 
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Fig. 3.12 Comparison of correlated and noncorrelated 
results in pure H 2 O for p- 0.0, L=0.1 m. 
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Chapter 4 

RADIATIVE INTERACTIONS IN LAMINAR 
FLOWS USING MONTE CARLO SIMULATION 

In order to further establish validity of the MCM, a relatively simple problem of 
radiative interactions is considered. The physical problem considered is that of steady- 
state energy transfer in laminar, incompressible, fully developed flow with constant 
properties in an absorbing-emitting gas between two parallel plates (Fig. 4.1). The 
condition of uniform surface heat flux is assumed such that the surface temperature 
varies in the axial direction. This problem is selected because gray as well as nongray 
solutions for this case are available in the literature [5, 65]. 


4.1 Basic Theoretical Formulation 


The energy equation for the presical physical system can be expressed as [8] 


dT dT\ _ d 2 T 
dH) _ dy* 




dp 

dx 


du \ 2 

dj) 


(4.1) 


where u and v denote the x and y components of velocity, respectively. In deriving 
Eq. (4.1), it has been assumed that the net conductive heat transfer and radiative heat 
transfer in the x direction (parallel to the plates) can be neglected in comparison to the 
flux variations in the y direction (normal to the plates). If, in addition, it is assumed that 
the Eckert number of the flow is small, then Eq. (4.1) reduces to 


dT dT 

4 - 7 ) 

dy 


P C p( u fc+ v —) = k 


d 2 T dg r 
dy 2 dy 


(4.2) 


The neglect of axial conduction and radiation in Eq. (4.2) is consistent with the 
formulation used in Ref. 9. 
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Fig. 4.1 Laminar flow between parallel plates with constant wall heat flux. 
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For a fully-developed flow, v=0, and u is given by the well-known parabolic profile 
as 

u ~ 6u m (£ - ( 2 ) ; ( = y/L (4.3) 

where u m represents the mean fluid velocity. Also, for the flow of a perfect gas with 
uniform heat flux, dT/dx is constant and is given by 

dT/dx = (2 aq w )/(u m Lk) (4.4) 


A combination of Eqs. (4.2)-(4.4), therefore, results in 


h d 2 T dq r 
dy 2 dy 




(4,5) 


Equation (4.5) is the governing energy equation for the parallel plates geometry. The 
boundary conditions for this problem can be expressed as 


T( 0) = T(L) = T w . —(y = L/2) = 0 (4.6) 

It should be noted that all boundary conditions given in Eq. (4.6) are not independent; 
any two convenient conditions can be used to obtain specific solutions. 

The radiative transfer term in the energy equation makes computation difficult because 
it turns the differential equation into an integro-differential equation. One exception is for 
the case of a gray medium. In this case, the equation for radiative transfer is expressed 
as [8] 


1 ^ = l« dTi 


(4.7) 


dy 2 4" rv " / 2" dy 

Equation (4.7) is a second order differential equation and, therefore, requires two bound- 
ary conditions. For black walls and T w j=T w2 , the boundary conditions for Eq. (4.7) 
become 

9r (I/2)=0; lq r (0) = Udq r /dy) 3=o 


(4.8) 
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In the present study, attention is directed on the MCM in solving the radiative transfer 
term for a gray as well as a nongray medium. Before going into a detailed numerical 
analysis of the energy equation including the radiative transfer term, it is essential to 
define the quantity of primary interest. 

For heat transfer in simple flow problems, the quantity of primary interest is the 
bulk temperature of the gas. For a fully-developed flow between parallel plates, this is 
expressed as 

1 

0 b = (T h - T w )/(q w L/k) = 6 f 0(0 (t - £ 2 K (4.9) 

0 

where q w — h(T w — Tf>), and h represent the equivalent heat transfer coefficient 
(W/cm 2 -K). 


4.2 Solution Procedure 


There are two levels to the numerical method proposed here. The first is concerned 
with the discretization and solution of the energy equation, while the second is due to the 
numerical evaluation of the radiative flux term that is included in the energy equation. 

The energy equation, Eq. (4.5), is discretized by a finite volume technique. The 
domain between two parallel plates is divided equally into N finite volume elements. For 
the ith finite volume elements, <!>Vi, a combination of Eqs. (4.5) and (3.1) results in the 
discretized energy equation as 


k 


Tj+i — 2T| + Tj — i 
Ay 


l2q w Ay 


(&-<?) 


-f -Qv-sVi + QA-sVi - QsVi = 0 


(4.10) 


where the conductive heat transfer is discretized by a central difference scheme and the 
radiative heat transfer consists of Q v ^ vi , Q A _£ Vi and Q^ V i te 11118 * The energy balance in 
each volume element results in a set of simultaneous equations equal to the total number 
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of finite volume elements. Each equation contains an unknown temperature which cannot 
be calculated independently and an iterative solution is necessary. 

Before solving the energy equation, the radiative energy interchange in each equation 
must be evaluated. In this study, the radiative terms Q v .£ vi , Qa-^vi Q^ Vi are 
simulated by the MCM. For the gray medium, the Monte Carlo formulations employed 
are from Refs. 6 and 7. For the nongray medium, the one-dimensional Monte Carlo 
formulations, as presented in Eqs. (3.7)-(3.9) and (3. 14), have been applied. 

In the Monte Carlo simulation, Q v £ Vi ^ Qa-Svi ^ obtained based on the assumed 
temperature distribution; a new temperature distribution can be calculated by solving the 
set of simultaneous equations given by Eq. (4.10). Two typical methods have been 
developed to solve the energy equation with the Monte Carlo simulation. In one of 
these methods [66, 67], convective and conductive heat transfer, as well as Q v _£ vi anc * 
Q a £ Vi , are calculated based on the assumed temperature distribution, a new temperature 
distribution is obtained from the term Q$ Vi in the energy equation. The numerical 
experiments conducted in this study indicate that this method has a high probability 
of producing divergent simulation and, therefore, it is not suitable for problems with 
large variations in optical length. In the other method [20], only Qv-Sv, “d Q A -<5vi 316 
calculated based on the assumed temperature distribution, a new temperature distribution 
which is included in the convective and conductive heat transfer terms, as well as Q £ Vi , ls 
obtained by solving a set of non-linear equations. This latter method was employed in the 
present study and the solution was obtained by using the NEQNF routine, which solves 
a system of non-linear equations in IMSL Library Package [68]. The change in local 
temperature in each iteration of the calculation is determined and when the maximum 
change is less than 10“^, the solution is considered to have converged. 

The radiative heat transfer can be calculated easily by the MCM, but the accuracy of 
the results obtained is affected by the number of the radiative energy bundles used in a 
calculation. If high accuracy is needed, it will be necessary to take longer computational 
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time even if a simple model is analyzed. Several methods are available to reduce the 
computing time and obtain higher accuracy. One of these methods, applied to the gray - 

gas, is the differential emissive power emission (DPE) method [66, 67]. In the DPE 
method, not only positive radiative energy bundles but also negative bundles are used, p 

and the number of energy bundles emitted from a gas element is proportional to the 

jfP : 

difference between emissive powers from two consecutive iterations. This treatment | 

does not change the physical processes of the Monte Carlo simulation. The proof of the 

K • 

equivalence of the DPE and regular methods is given in the cited references. K - 

l 

4.3 Results and Discussion 

Based on the theoretical and numerical analyses described in the previous sections, 
a computer code, which is given in Appendix B, has been developed to investigate gray 
as well as nongray radiation interactions in incompressible flows between two parallel 
plates. For the case of black walls, gray analytical solutions and nongray approximate 
solutions (based on the method of variation of parameters) are available in the literature 
[5, 65]. In this study, the Monte Carlo solutions have been compared with these results 
for identical conditions. The absorbing-emitting media considered were pure H 2 O and > 

CO 2 . The results are expressed in terms of the non-dimensional bulk temperature. The 
plate spacings considered range from 0.01 cm to 100.0 cm. The calculation was carried : 

out on a Sun Workstation. The domain was divided into 40 finite volume elements with 
equal thicknesses. The total number of energy bundles selected was 50,000 for nongray ; 

and 200,000 for gray simulations. The amount of energy per bundle depends on the 
temperature. One of the important parameters related to the temperature distribution 
is the heat flux from the plates; care should be taken to choose this heat flux. In the 
solutions from the literature [5, 65], the assumption of linearized radiation was made and 

V, : { 

the radiative properties were considered to be independent of temperature. In order to ^ 

facilitate comparisons between the Monte Carlo solution and the approximate solution, 

£. : 
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different values of heat flux at the wall were chosen when the plate spacings were changed. 
The CPU time requirement for a converged solution with a specific plate spacing was 
on the order of ten seconds for the gray case if the DPE method was applied and on the 
order 1000 seconds for the nongray case. The numerical experiments conducted in this 
study indicate that the DPE method can reduce the CPU time about an order of magnitude 
compared to the regular method without loss in the accuracy of results. 

Figures 4.2-4.5 show comparisons between the gray analytical solutions and the gray 
Monte Carlo solutions for different media, temperatures, and pressures. The medium 
considered is C0 2 in Figs. 4.2 and 4.3. In Fig. 4.2, the pressure of C0 2 was kept 
at 1.0 atm but the plate temperatures were 500 and 1000 K. In Fig. 4.3, the wall 
temperature was kept at 1000 K but the pressures was changed from 1.0 to 5.0 atm. 
The figures show that the predictions by the MCM are very close to the analytical 
solutions at different temperatures and pressures. Figures 4.4 and 4.5 show the results 
for H 2 0. Similar to the case for C0 2 , the Monte Carlo solutions were found to be in 
good agreement with the analytical solutions in the H 2 0 medium at different temperatures 
and pressures. The results demonstrate that radiative interactions are enhanced and the 
temperature distribution becomes more uniform between the parallel plates with increases 
in temperature and pressure. 

Figures 4.6-4.9 show comparisons between the nongray approximate solutions based 
on the method of variation of parameters and the nongray Monte Carlo solutions for 
different media, temperatures, and pressures. The medium considered is C0 2 for the 
results presented in Figs. 4.6 and 4.7. In Fig. 4.6, the pressure of C0 2 was kept at 
1.0 atm but plate temperature was changed from 500 to 1000 K. In Fig. 4.7, the wall 
temperature was kept at 1000 K but the pressure was varied from 1.0 to 5.0 atm. The 
figures show that the Monte Carlo solutions compare favorably with the approximate 
solutions at different temperatures and pressures. Figures 4.8 and 4.9 show the results 
for H 2 0. Similar to the case of C0 2 , the Monte Carlo solutions essentially match the 
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Fig. 4.3 Comparison of gray solutions for CO 2 at T w =1000 K. 











Fig. 4.6 Comparison of nongray solutions for CO 2 at p=l atm 
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analytical solutions for H2O at different temperatures and pressures. The effects of an 
increase in temperature and pressure on the radiative interactions and the temperature 
distributions between the parallel plates in the nongray cases are also found to be similar 
to those in the gray cases. 

The results presented in this chapter demonstrate clearly that the one-dimensional 
nongray Monte Carlo formulatons developed in the previous chapter are very reliable 
and accurate. These formulatons have been also applied to investigate the radiative 
interactions in entry region turbulent flows, and detailed information is available in Ref. 
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Chapter 5 

MONTE CARLO SIMULATION FOR RADIATIVE 
TRANSFER IN MULTI-DIMENSIONAL SYSTEMS 

In Chap. 3, radiative heat transfer between two infinite parallel plates was simulated 
in an exact manner. However, application of this exact treatment to multi-dimensional 
problems can be extremely complicated and numerical solutions to these formulations can 
be very difficult. However, by introducing an appropriate assumption, the complicated 
Monte Carlo formulations in multi-dimensional problems can be simplified significantly. 
In this chapter, attention is directed to a two-dimensional problem. The physical problem 
is described in Sec. 5.1. The exact Monte Carlo formulations are developed in Sec. 
5.2. The approximate Monte Carlo formulations are developed in Sec. 5.3. Comparisons 
between exact and approximate Monte Carlo solutions are made in Sec. 5.4. 

5.1 Physical Problem 

Consider an absorbing and emitting molecular gas between two parallel plates of 
finite length L and height H and infinite width, as shown in Fig. 5.1. The inlet and outlet 
of the gas are at x=0 and x=L, respectively, and both ends are treated as pseudoblack walls 
with prescribed temperatures. Temperature, concentration and pressure in the medium are 
assumed to be known. The walls are assumed to be diffuse but not necessarily gray. The 
wall temperature distribution is also specified. In order to calculate the radiative source 
term —V.q r inside the medium and the net radiative wall flux q rw , the medium considered 
is divided into an MXxMY array of rectangular volume elements (Fig. 5.1). Similarly, 
the two real walls are each divided into MX surface elements, and the inlet and outlet 
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pesudo walls are each divided into MY surface elements. Temperature, concentration 
and pressure are assumed to be constant in each element. 

The following Monte Carlo analysis is based on an arbitrarily chosen finite volume 
element ABCD (Fig. 5.2) with the length and height equal to b and c, respectively. Exact 
correlated and non-correlated formulations are derived first; then approximate correlated 
and non-correlated formulations are developed. The statistical relationships for an energy 
bundle emitted from a surface element in each case can be derived by following the 
same procedure. 

5.2. Exact Correlated and Non-correlated Monte Carlo Formulations 

In this case, an energy bundle is simulated in an exact manner in terms of the 
narrow band model without approximation. Let us consider the Planck spectral blackbody 
intensity Ibo> that enters the element ABCD at some point s on side AB and intersects 
one of the other three sides of the element at the point s', as shown in Fig. 5.2. It should 
be understood that each side of the element is a surface. A spherical coordinate system 
is established and centered at the point s. The distance between the points s and A is 
x*. From Ref. 6, the amount of energy emitted in the wavenumber interval duo, along a 
pencil of column s-~>s' with a solid angle increment dfl and an area increment dx* is 

dQ = hw [l - t w (s — ► s')] cos 9d£ldx*du> (5.1) 


The total emitted energy, calculated in terms of the intensity entering from the sides of AB 
(Oc0<tt) and DC (7r<0<27r), is obtained by integrating Eq. (5.1) over the wavenumber. 


polar angle, azimuthal angle and area as 

oo b ic 2ir 

o-JJJJ Ibuj [l — t w (s —> s)] cos 9 sin OdifrdOdx* duo 

oooo 

Referring to Fig. 5.2, the distance ss' is expressed as 

min {cf cos 6, (b — rr*)/(cos7/>sin0) }, — 7 t/2 < < 7t/2 

ss' — 


(5.2) 


min{c/ cos 9 , 


x*/(cosijjsm9)}, / K/2 <\j> < Stt/2 


(5.3) 





63 


The value of ss' cannot be calculated from just one expression because the point s' may 
be located on different sides of the element ABCD. All the possible travelling paths of 
the intensity in the element ABCD should be considered to evaluate the value of Q. 

Similar procedures can be used to obtain expressions for the emitted radiative energy 
calculated in terms of the intensity entering from the sides of AD and BC. Thus, the total 
emitted radiative energy from the finite volume element ABCD consists of two terms. 
They represent the emitted energy calculated in terms of the intensities entering from sides 
AB, DC and from sides AD, BC, respectively, and cannot be manipulated algebraically 
into one term. Usually, the statistical relationships for simulating an energy bundle 
emitted from a volume element in the MCM are developed from the formulation for the 
total emitted radiative energy from this volume element. However, this can complicate 
the analysis since there exists two independent terms in the formulation of total emitted 
radiative energy. In this study, the two independent terms are treated separately, and 
the Monte Carlo analysis is based on a single term. This means that the Monte Carlo 
analysis is based on Eq. (5.2) if an energy bundle in element ABCD starts from either 
side AB or DC. Otherwise, the Monte Carlo analysis is from another term. 

The Monte Carlo formulations presented here are developed on the basis of Eq. 
(5.2). The simulation of an energy bundle includes the determination of wavenumber, and 
starting point and direction of emission of this energy bundle in the finite volume element. 
The statistical relationships for determining these parameters are obtained readily from 
Eq. (5.2) as [6, 54, 55] 


R 


u) 


b ir 2 tt 

III I hull- 


0 0 0 0 


r w (s 


— > s 1 )] cos 9 sin 6dtf>d0dx* doj 

~Q 


(5.4) 
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Re 
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(5.6) 
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ip oo b ic 
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0 0 0 0 


— ► s')] cos 9 sin 9d9dx* dwd\p 

~Q 


(5.7) 


In Eqs. (5.2) and (5.4)-(5.7), r w is a real spectral transmittance. Before solving these 
equations to obtain u>, x\ 9 and x/) from a set of given values of Ru,R x *,R0,Ril>> the 
narrow band model should be applied to approximate the real spectral transmittance. 

Taking the spectral average over all narrow bands and using the narrow band 
approximation as in Eq. (3.10), Eqs. (5.2) and (5.4)-(5.7), are expressed as 


( b ic 2 t 
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( 5 . 12 ) 


Similar to the one-dimensional problem analyzed in Chap. 3, in order to solve Eqs. (5.8)- 
(5.12) for a set of given values of Ra;, R^, R# and R^, interpolation and approximation 
methods must be employed. 
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Knowing wavenumber, emission point and travelling direction of an energy bundle, 
the next question is where the energy bundle is absorbed. Let us still consider the emitted 
radiant energy along a pencil of column s— >s' (Fig. 5.1). After this amount of energy is 
transmitted over a column s'— >s", the remaining radiant energy is given by 


dQ f = Ifa [l — Ta,(s — > s')] 7 -^( 5 ' — > s") cos 9dCldx*dcj 


(5.13) 


where r w (s'— ►s") is the spectral transmittance over the path s'— >-s". Taking a narrow 
band average of Eqs. (5.1) and (5.13) and dividing the latter by the first, the statistical 
relationship for determining the location of absorption can be expressed as 


Rl = 


[1 - t u (s -> -» s") 

1 - %;{s s') 


T^s' —> s") ~ T w (s — ► S , ) r w( s/ — > s") 
1 - -► ^') 


(5.14) 


Similar to Chap. 3, the averaged product t u (s — > s , )r a ,(s' — > s") can be treated in a 
spectrally correlated or non-correlated manner. The first choice results in the spectrally 
correlated formulation as 


Ri = 


Tu>(s’ -+ s”) - ^(5 -» s") 

1 - 7^(S -> S') 


(5.15) 


and the latter choice results in the spectrally non-correlated formulation as 


Rl = - s") 


(5.16) 


Therefore, it is seen that exact correlated and non-correlated Monte Carlo formulations 
differ only in the relation for Ri as given in Eqs. (5.15) and (5.16). 

Comparing each formulation in this chapter with the corresponding one-dimensinal 
formulations developed in Chap. 3, it is found that the exact correlated and non-correlated 
statistical relationships for Ri are the same but statistical relationships for R u , R x * , Re, R^ 
are different. This phenomenon is true for any two different problems. 
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53 Approximate Correlated and Non-correlated Monte Carlo Formulations 


In the exact analysis, the Monte Carlo formulations need to be developed from each 
independent term in the expressions for the total emitted radiant energy of a volume 
element. Numerical evaluations of Eqs. (5.8M5.12) for Q, Rl>-> Rx*, R$, involve 
four-dimensional integrations and the integrands in these equations ar complex functions 
of integration variables. Obviously, the Monte Carlo simulation is already complicated 
although the problem considered is a simple two-dimensional problem. The difficulty may 
continue to increase considerably if the complexity of the problem increases. To simplify 
complicated Monte Carlo formulations, it is assumed that the volume dV of a volume 
element is very small so that the energy emitted within dV escapes before reabsorption. 
This assumption has been used widely in many studies to simplify radiation analysis. 
The total emitted radiative energy and the statistical relationships for determining the 
wavenumber and emission direction of an energy bundle from a finite volume dV are 
given by [6, 54, 55] 

oo 

QdV = 4tt J hudV du) (5.17) 

o 

w 

Ru, = 75 < 5 * 18 ) 

f 
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Re = 


1 — cos 0 
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(5.19) 


R* = £ ( 5 - 2 °) 

The emission point of an energy bundle from a volume element is assumed to be the 
center point of the element. This assumption is justifiable in an infinitesimal volume 
element. Introducing the narrow band approximation, Eqs. (5.17) and (5.18) become 

m w 

QdV = ^ (l^I bu ,kAuj k ^dV 

k = 1 


(5.21) 
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4?r ^2 (K w kI bu; kAuj k )dV 
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k = i 


Qdv 


(cj n_1 <u<u n ) 


(5.22) 


Here, k u k is the mean absorption coefficient over a narrow band and is obtained as [29] 

„ _ In T„*(Lm) (5.23) 

L m 

where Lm is the mean beam length of the volume element. It is evident that Eqs. (5.19)- 
(5.22) are much simpler than the corresponding equations for the exact treatment of the 
Monte Carlo simulation. These simple formulations do not change with the complexity 
of the problem. 

The statistical relationship presented here for determining the location of absorption 
of an energy bundle emitted from the volume element dV is different from that available 
in literature [6, 54, 55]. This is because a narrow band model is incorporated in the 
present formulation. Equation (5.15) is the general formulation to calculate Ri with 
consideration of the spectral correlation. Substituting the mean transmittance with the 

mean absorption coefficients in the denominator of Eq. (5.15), yields 

TZ(s' -► s") - s") 
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Since dV is very small, the following approximation can be invoked 
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This approximation is also applied in deriving Eqs. (5.17)-(5.20). Consequently, Eq. 
(5.24) is simplified as 
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Equation (5.26) is the approximate correlated statistical relationship for determining 
the location of absorption and it is different from the corresponding exact correlated 
formulation as given in Eq. (5.15). The approximate non-correlated statistical relationship 
for determining the location of absorption cannot be simplified further and it is the same 
as given by Eq. (5.16). Therefore, similar to the exact correlated and non-correlated 
formulations, the approximate correlated and non-correlated formulations differ only in 
the expression for R\. 


5.4 Results and Discussion 


In order to validate the approximate Monte Carlo analyses and to investigate the 
effects of spectral correlation, two problems have been selected by referring to the work 
of Zhang et al. [28]. The results for the net radiative wall flux and the radiative source 
term have been obtained for four different formulations which correspond to the exact 
correlated solution, approximate correlated solution, exact non-correlated solution and 
approximate non-correlated solution. In the problems considered, the length and height 
of two parallel plates are L=1.2 m and H=0.6 m, respectively. The two wall emissivities 
are chosen to be the same and equal to 0.8. The total pressure of the gas is taken to 
be 1 atm. One of the problems considered is an isothermal and homogeneous H 2 O-N 2 
mixture in which the mole fractions are: X# 2 o=0.6 and Xn 2 =0.4\ the gas temperature is 
1500 K; and the real and pesudo walls are held at 300 K. The other problem considered 
is a nonisothermal and inhomogeneous H 2 O-O 2 -N 2 mixture in which the mole fraction 
distributions are given by 

Xh 2 o = 0-3 

*o 2 = 0.1 

Xn 2 = 1 - XR 2 0 - x 02 

and the gas temperature distribution is assumed to be 

T(x,y) = 1000 + 1200 
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The two real walls and the inlet pseudo wall are kept at a temperature of 1000 K. The 
outlet of the gas is open to a 300 K atmosphere, so the temperature of the outlet pseudo 
wall is 300 K. For both problems, only H 2 O is considered to be a radiatively participating 
species. There are five important absorption bands for H 2 O. All of these bands have been 
taken into account in this study and they consist of mo?=295 narrow bands in the range 
from 150 cm -1 to 7500 cnr 1 . 

To assure that the statistical results make sense in the Monte Carlo simulation, two 
requirements must be met. One is the accuracy of the statistical results for a given 
grid. The other is the independence of the results from the grid. In this study, the 
designated statistical accuracy of the results is defined in such a way that when the 
relative statistical errors are less than ±5%, the probability of the results lying within 
these limits is greater than 95%. Independence of the results on a grid is considered 
to have been achieved when the medium is divided into 20x20 uniform finite volume 
elements for the problems considered. For this grid, the total number of energy bundles 
had to be 2,000,000 in order to meet the designated statistical accuracy requirement. All 
calculations have been carried out on a Sun Sparc Workstation. The CPU times required 
for different solutions for two different problems are listed in Table 5.1. It should be 
noted that the present computer code was written for problems involving nonisothermal 
and inhomogeneous mixtures. No efforts have been made to simplify the problem for an 
isothermal and homogeneous mixture specifically. For integrations and interpolations 
in the exact Monte Carlo formulations, Eqs. (5.8)-(5. 12), the divisions of the side 
length, polar angle and azimuthal angle (within half of their ranges in a rectangular 
volume element) were chosen to be m x *=10, m^=10 and m^,=10, respectively. The 
emitted radiative energy from each of the mo;xm x *xm^xm^=295x 10x10 x 10 medium 
columns was then calculated and stored. The required integrations and interpolations were 
implemented from the summation of the values of radiative energy in different columns. 
These computations were done for each volume element. Obviously, this procedure is 
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Table 5.1 CPU time (minutes) required for different solutions 



Exact 

correlated 

Approximate 

correlated 

Exact 

non-correlated 

Approximate 

non-correlated 

Isothermal and 

homogeneous 

mixture 

265 

112 

325 

170 

Nonisotherm al 
and 

inhomogeneous 

mixture 

269 

167 

378 

225 
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very time consuming. This is the major reason why the CPU time for exact solutions 
is much larger than that for approximate solutions in Table 5.1. It should also be noted 
that determination of the absorption location of an energy bundle, by using Eq. (5.15) 
and Eq. (5.26), takes about the same amount of time. 

The problem with an isothermal and homogeneous mixture is considered first. The 
behavior of four different solutions is illustrated in Figs. 5.3-5.5. Figures 5.3 and 5.4 
show radiative source distributions at locations equal to x/L=0.225 and x/L=0.5 on the 
plates, respectively. The approximate correlated results agree with the exact correlated 
results and the approximate non-correlated results agree with the exact non-correlated 
results. As the distance from the walls increases, all four solutions predict the same trend 
in the radiative source results. The two non-correlated solutions are far below the two 
correlated solutions. 

The distribution of radiative wall heat flux along the plates is presented in Fig. 5.5. 
The approximate correlated solution is found to be almost the same as the exact correlated 
solution and the approximate non-correlated solution is seen to be slightly higher than the 
exact non-correlated solution. The difference between the correlated and non-correlated 
results is seen to be significant. For the most part, the two non-correlated solutions are 
approximately two times higher than the two correlated solutions. 

The results for a nonisothermal and inhomogeneous mixture are illustrated in Figs. 
5.6-5.9. The H 2 O mole fraction calculated from Eq. (5.27) has a maximum value at 
the mid plane of the geometry considered, and decreases gradually away from the center 
point The temperature in the medium, calculated from Eq. (5.28), increases away 
from the walls and the inlet. Figures 5.6-5.8 show the radiative source distributions 
at locations equal to x/L=0.275, 0.5, and 0.825 along the plates, respectively. As the 
distance from the inlet location increases, the temperature change becomes more steep 
and temperatures in the central region are high. Thus, the change in radiative source 
results is becoming abruptly as seen from Figs. 5.6-5.8. In all three figures, it is 








Fig. 5.4 Radiative source distribution at the location 
x/L=0.5 for isothermal and homogeneous H 2 O-N 2 mixture. 
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Fig. 5.6 Radiative source distribution at the location x/L= 0.275 
for nonisothermal and inhomogeneous H2O-O2-N0 mixture. 
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Fig. 5.7 Radiative source distribution at the location x/L=0.5 
for nonisothermal and inhomogeneous H 2 O-O 2 -N 2 mixture. 
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Fig. 5.8 Radiative source distribution at the location x/L=0.825 
for nonisothermal and inhomogeneous H 2 O-O 2 -N 2 mixture. 
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Fig. 5.9 Radiative wall flux distribution for 
nonisothermal and inhomogeneous H2O-O2-N2 mixture. 
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evident that the approximate correlated solution is in good agreement with the exact 
correlated solution and the approximate non-correlated solution approximates the exact 
non-correlated solution. The approximate correlated solution appears to be slightly higher 
in the wall region and lower in the central region than the exact correlated solution. The 
difference between the correlated and non-correlated solutions is significant as in the first 
problem. From the correlated solutions, it is evident that the gas goes from a net absorber 
near the walls to a net emitter away from the walls. On the other hand, the non-correlated 
solutions predict that the gas is a net emitter in nearly all regions. 

Figure 5.9 illustrates the distribution of radiative wall flux along the plates. The 
radiative wall flux is seen to increase at first, reach a peak value near the outlet, and then 
decrease. Such behavior is due to the fact that, for the problem considered, the outlet 
region is equivalent to a cold sink. This cold sink has a strong effect on the radiative 
heat transfer in the adjacent region. Among the four different solutions, the approximate 
correlated solution is slightly lower than the exact correlated solution and the approximate 
non-correlated solution is slightly higher than the exact non-correlated solution. A 
comparison of different solutions reveals that the non-correlated formulations predict 
much higher radiative energy absorption on the walls than the correlated formulations. 
The difference in results can be reach as high as one order of magnitude at some locations. 

From the results presented, it is evident that approximate formulations can provide 
results very close to those from the corresponding exact formulations. The non-correlated 
formulations, however, predict much lower radiative source distributions in the medium 
and much higher radiative wall fluxes along the plates than the correlated formulations. 
The reason for this difference is the same as that for one-dimensional problem. That is, the 
Ri calculated from the non-correlated formulation, Eq. (5.16), is greater than that from 
the correlated formulations, Eqs. (5.15) and (5.26). Therefore, for the non-correlated 
formulation, an energy bundle travels a long distance and is likely to be absorbed on the 
wall. This also explains why the CPU time required for the non-correlated solution is 
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larger than that required for the corresponding correlated solution (Table 5.1). Because 
of significant differences between the correlated and non-correlated solutions, the same 
conclusion as that in Chap. 3 is drawn that the non-correlated formulations are not useful. 
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Chapter 6 

RADIATIVE INTERACTIONS IN CHEMICALLY 
REACTING COMPRESSIBLE FLOWS 

In Chap. 5, an accurate radiation transport model using the approximate Monte Carlo 
correlated formulations, has been developed and validated. The formulations in this model 
are simple and can be applied easily in nongray and multi-dimensional systems. The 
objective of this chapter is to apply the approximate Monte Carlo correlated formulations 
to investigate the radiative interactions in multi-dimensional chemically reacting flows. 
The basic formulations are provided in Sec. 6.1. The method of solution is presented in 
Sec. 6.2, and the results and discussion are contained in Sec. 6.3. 

6.1 Basic Formulations 


6.1.1 Physical Model 

As mentioned in the introduction, there has been extensive research directed toward 
the development of scramjet propulsion systems. To investigate the radiative effects on 
these systems, a specific physical model will be considered in this study which is a 
supersonic flow of premixed hydrogen and air in an expanding nozzle (Fig. 6.1). The 
nozzle wall is modeled, as noted, by a shifted sinusoidal curve. The inlet temperatures of 
hydrogen and air are considerably high so that chemical reactions take place in the entire 
flowfield. The products of hydrogen-air combustion include water vapor and hydroxyl 
radicals. These species are highly absorbing and emitting. To simulate the flowfield 
accurately, all important phenomena such as chemistry, radiation and turbulence should 
be taken into account and the fully elliptic form of the governing equations must be used. 
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6.1.2 Governing Equations 


In this study, the two-dimensional nozzle flow considered is described by the Navier- 
Stokes and species continuity equations which can be represented in the physical coor- 
dinates as 

3U OF dG „ 

ir + ^ + ^ = H (6A) 

where vectors U, F, G and H are given by 


U = 


pu 

pv 

pE 

Lp/i-l 


pu 


G = 


pu — (T x 

puv — Tyx 

(, pE - 0 x )u — T xy v + q x 

pfi(u + Ui) 
pv 

pUV - T xy 
pv 2 - cry 

(pE — <Ty)v - Tyx U + q y 
pfi(v + Vi) 

-o 
0 
0 

-Wi 


V*9r 


( 6 . 2 ) 


(6.3) 


(6.4) 


(6.5) 


84 


The other terms appearing in vectors F, G, and H are defined as 

, [du dv\ du 

^ = - p+X [d- X + Yy) +2tl Yx 

, ( du dv\ dv 

^ = ~ p + A (^ + %) +2/t % 


( 6 . 6 ) 


(6.7) 


/ du dv 

Txy = Tyx= V\fa + fy 


dT 


N, 


q x 


q y 


= -k-^+pY, h <f& 


dT 


*=i 


n 3 


= - k o^+pY, h 'f<V' 


x'=i 


2,2 

E = .l + 1^L + Ehifi 

9 L i=l 


( 6 . 8 ) 


(6.9) 


( 6 . 10 ) 


( 6 . 11 ) 


k = h? + J 


C Pi dT 


t r 


N a 


P 


i=l 


■*'Ei 


where A = —\p. 


( 6 . 12 ) 


(6.13) 


In Eqs. (6.1), only (N s — 1) species equations need to be considered since the mass 
fraction of the species is prescribed by satisfying the constraint equation 


N a 

£/. = > 


(6.14) 


t=l 


The diffusion velocity of the ith species is obtained by solving the Stefan-Maxwell 
equation [70], neglecting the body force and thermal diffusion effects, as 


VX, 




(6.15) 
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The preceding equation is also applied only to (N s — 1) species. The diffusion velocity 

N a „ 

for the remaining species is prescribed by satisfying the constraint equation fiVi = 0, 

t=i 

which ensures the consistency. 

6.1.3 Thermodynamic Model 

To calculate the required thermodynamic quantities, the specific heat for each species 
C Pi is first defined by a fourth-order polynomial in temperature, 

% = .4; + B,T + CiT 2 + AT 3 + EiT 4 (6. 16) 

R 

The values of the coefficients appearing in the equation are found in Ref. 71. Knowing 
the specific heat of each species, the enthalpy of each species can be found from Eq. 
(6.12) and the total internal energy is computed from Eq. (6.11). 

6.1.4 Chemistry Model 

Chemical reaction rate expressions are usually determined by summing the contribu- 
tions from each relevant reaction path to obtain the total rate of change of each species. 
Each path is governed by a law of mass action expression in which the rate constants 
can be determined from a temperature dependent Arrehenius expression. In vector H, 
the term w{ = M{C t represents the net rate of production of species i in all chemical 
reactions and is modelled as: 

E^rE^ ; j = h---N r (6.i7) 

t=i bj t=i 

N r [ N, N, 

W, = M,c, = MiJ2 ( 7 ." - 7 ,';) kf, II <#” - k h II C ™ m < 6 - 18) 

j—\ L m=l m—1 

Equation (6.17) represents an N r step chemical reaction and Eq. (6.18) is the production 
rate for the ith species. The reaction constants kf } and kf )] are calculated from the 
following equations: 

k fl = AfT^expf-^L) ; j = l,---N r (6.19) 



86 


kbj kfj/keqji j 1? * ' * N r 


( 6 . 20 ) 


The equilibrium constants appearing in Eq. (6.20) are given by 

_ (j_ 


where 


keq ’ 1 R U T ) eXP { R U T 


j = !,■■■ N r 


N, 


N, 


A = j = h---Nr 


1=1 1=1 


( 6 . 21 ) 


( 6 . 22 ) 


N a Ns 

A G Rj = 'fijSi - Y 't'! 9 ’’ j = h-" Nr 

t=l t = l 


(6.23) 


U- = AAT - InT) - ~T 2 - ^T 3 - ^T 4 
Ri v ’ 2 6 12 


-|j-T 5 + F { - GiT ; 


i = - ■ N r 


(6.24) 


The forward rate for each reaction is determined from Eq. (6.19). The hydrogen- 
air combustion mechanism used in this work is from Ref. 3, but only seven species 
and seven reactions were selected for this study. The constants Aj, Nj and Ej for these 
reactions are listed in Table 6.1. The species Gibb’s free energy expression Eq. (6.24) is 
obtained from the integrations of the specific heat C Pi and the coefficients in Eq. (6.24) 
are obtained in the same way as in Eq. (6.16). 

6.1*5 Diffusion Models 

The viscosity , thermal conductivity, and diffusion coefficient consist of the contri- 
butions from both fluid molecules and turbulent flow and they are expressed as 


fi = /q + 
k = k\ + k t 
Dij = D\j + D\ 3 


(6.25) 
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Table 6.1 Hydrogen-Air Combustion Mechanism (7 species, 7 reactions) 


No. 

Reaction 

A 

N 

E 

m 

H 2 + 0 2 -► OH + OH 

1.70E+13 

0.0 

24233 

2 

H + 0 2 -» OH + 0 

1.42E+14 

0.0 

8250 

3 

OH + H 2 -4 H 2 0 + H 

3.16E+07 

1.8 

1525 

4 

0 + H 2 -> OH + H 

2.07E+14 

0.0 

6920 

5 

OH + OH — > H 2 0 + 0 

5.50E+13 

0.0 

3523 

6 

H + OH + M -4 H 2 0 + M 

2.21E+22 

- 2.0 

0 

7 

H + H 4 - M — > H 2 + M 

6.53E+17 

- 1.0 

0 
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where ///, h, D\- represent the molecular mixture viscosity, thermal conductivity, and 
diffusion coefficient, respectively; fi t , k t , D\- represent the turbulent viscosity, thermal 
conductivity, and diffusion coefficient, respectively. 

The individual species molecular viscosities are computed from Sutherland’s law 


jh _ = ( t y V2 ?!„ + .<?■ 
f*0i \Tq i ) T + Si 


(6.26) 


where //oi and Toi are reference values and So is the Sutherland constant. All three 
values are tabulated for the species in Refs. 72 and 73. Once the molecular viscosity 
of each species has been determine, the molecular mixture viscosity is determined from 
Wilke’s law [74] 




N a 

H = £ 

*= 1 1 + 37 X) Xjfaj 

j= iUM) 


N a 


(6.27) 


where 


fiij — 


{l + [(w/w)(«/«)] I/2 (W i /J#i) 1/4 }' 

4y/2[l + (Mi/Mj )] 1 / 2 


(6.28) 


The individual species thermal conductivities are also computed from Sutherland’s 


law 


h 

hi 




t \ 3/2 r', + s[ 


(6.29) 


T'J T + S'i 

but with different values of the reference values hi and T oi and the Sutherland’s constant 
Sj. These values are also taken from Refs. 72 and 73. The molecular mixture 
thermal conductivity is computed using conductivity values for the individual species 
and Wassilewa’s formula [75], 

\r 

h 


N s 

*i = £ 


(6.30) 


1-1 1 + x X) 

j= iUM) 

where 4 = 1.075^ and <f>ij is taken from Eq. (6.28). 
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For dilute gases, Chapman and Cowling [70] used kinetic theory to derive the 
following expression for the molecular binary diffusion coefficient Djj between species 
i and j, 

D , = 0.001858T 3 / 2 [(Mi + ] 1/2 (6 3J) 

tj pvjjttD 

Here, the diffusion collision integral ft# is approximated by 

Cl D = r*-°- 145 + (T* + 0.5) -2 (6.32) 


where T* = T/T £ij . The values of the effective temperature Tsij and effective collision 
diameter cr^ are taken to be averages of the separate molecular properties of each species, 
giving [70] 

a ij = + a j) < 6 - 33 ) 


and 



(6.34) 


To evaluate the turbulent viscosity pt, a turbulence model needs to be selected. An 
appropriate model selected in this study is the Baldwin-Lomax model. This model is 
very convenient to use and is also reliable for the flows like those considered here. The 
description of this model can be readily found in the literature [76-78]. Knowing turbulent 
viscosity pt, the turbulent thermal conductivity kt and turbulent diffusion coefficient D\- 
are calculated from the turbulent Prandtl number and the turbulent Schmidt number, 
respectively. 

6.1.6 Radiative Transfer Model 

The radiative effects on the nozzle flowfield arise through the term —V.q r in the 
energy equation and the radiative effects on the heat transfer on the nozzle walls arise 
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through the term q rw • The exact expressions for both —V.q r and q rw are very complicated 
integro-differential equations and they are usually treated separately from the governing 
equations. Therefore, the term —V.q r has been moved to the right hand side to be taken 
as the source term in the energy equation. As indicated earlier, the approximate Monte 
Carlo correlated formulations as seen in Eqs. (5.19)-(5.22) and (5.26) are employed 
to simulate the radiative heat transfer term. This treatment can provide a quantitative 
prediction of radiative interactions for the present problem. 

6.2 Method of Solution 


6.2.1 Grid Generation 

Equation (6.1) is written in the physical domain (x, y) and must be transformed to 
an appropriate computational domain (£, p) for solution. An algebrabic grid generation 
technique developed by Smith and Weigel [79] was used for grid generation in this study. 
From the computational point of view, it is desirable to have a uniform rectangular grid 
enclosed in a cube, where the exterior of the cube represents the physical boundaries. 
To have such grids, a body-fitted coordinate system was transformed linearly from the 
physical domain (x, y) to the computational domain (f, p) as follows: 

£i = x(£,0) Lower 

yi=y(p,0) Boundary (6.35) 


®2 = s(fjl) Upper 

y 2 =y(p,l) Boundary (6.36) 


x = x(£,l)p + x(£, 0)(1 - 77 ) Between the 

y = y({>i)>7 + y({»o)(i -q) 


Boundaries 


(6.37) 
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where 0 < f < 1; 0 < // < 1. The grid should be concentrated in the regions 

of high gradients to predict the solutions accurately. Therefore, more grid points are 
required near the solid boundaries. The concentration of the grid in the //-direction can 
be accomplished by 

- _ {0y + 1) ~ (P* ~ 1) exp [~C(i? - 1 + <*)/( 1 - a)] 

^ (2a + 1){1 + exp[C(»/ - 1 + a)/(l - a)]} 

where 

c -‘"(s ± r) «’> 

If a is equal to zero (c*=0), the compression takes place only near the lower wall ( 77 = 0), 
and if a is equal to one half (< 2 = 1/2), the compression takes place near both walls. The 
term {3 y has a value between one to two, and as it gets closer to one, the grid becomes 
more concentrated near the walls. Employing this concentration, Eq. (6.37) is written 
in terms of fj as 

* ' = + a?(f,0)(l -fj) 

y = y(£, 1)7 + y(f , o)(i - fj) (6.40) 

where 0 < rj < 1 . 

Based on the above analysis, the grid mesh for the present problem is generated as 
seen in Fig. 6.2. Because the flow is assumed to be symmetric about the centerline of a 
two-dimensional nozzle, only the upper half of the nozzle is shown. It should be noted 
that the grid is concentrated in the normal direction in order to capture the boundary 
layer and the grid is kept uniform in the flow direction. 

The above grid mesh was used for the flowfield simulation, but, the grid mesh for 
radiation simulation was quite different. A uniform grid mesh as seen in Fig. 6.3 was 
applied for radiation simulation for the present problem. Such a grid mesh is justifiable 
because radiation is a long-range phenomena and there is no need to use a concentrated 
grid mesh. 
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Fig. 6.3 Grid mesh for radiation simulation. 
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6.2.2 Numerical Algorithm 


The governing equations, Eqs. (6.1), are expressed in the computational domain as 

dU dF dG - 

H + l6 ' 41) 

where 


U = UJ 
F = Fy v - Gx v 
G = Gx{ ~ Fy^ 
H = HJ 


J = - V£X v 


(6.42) 


Here x^, x v , y^, y v are the transformation matrices and J is the Jacobian of the transfor- 
mation. The matrices can be computed numerically once the physical grid coordinates 
have been prescribed. 

The governing equation system, Eq. (6.41), can be stiff due to the kinetic source 
terms contained in the vector H. To deal with the stiff system, the approach used in Refs. 
80 and 81 was followed and the kinetic source terms were computed implicitly. In a 
temporally discrete form, Eqs. (6.41) then become 


{/ n+ 1 = U n -At 


dF 

dt 


+ 


dG 

drj 




(6.43) 


After employing a Newton linearization for H and rewriting in delta form, Eq. (6.43) 
becomes 


[/ - MK n )&U n + l = -A tR n 


(6.44) 


where 


dF 

* n ='ar 


+ 1^1 -H n 
dr) 


(6.45) 


is the steady-state residual, I is the identity matrix, K n is the Jacobian of H with respect 
to U, ( dH/dU ), and A U n+1 = U n+1 - U n . 
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Once the temporal discretization used to construct Eq. (6.44) has been performed, 
the resulting system is spatially differenced using the explicit, MacCormack predictor- 
corrector schemes [82]. This results in a spatially and temporally discrete, simultaneous 
system of equations at each grid point [80]. Each simultaneous system is solved using 
the Householder technique [81] in combination with the MacCormack technique, which 
is then used to advance the equations in time. The modified MacCormack technique 
then becomes 


[/ - Af/<7‘]A?7" + ' = -AtS + R”j 
t/jp = Ufj + A t/£+* 

I - Af/qp] A//" +1 = - AM'/?" 
i/£ +1 =U n + 0.5 [A //£+* + At/” +1 ] (6.46) 

where S + R represents a forward spatial difference of R and S~R a backward spatial 

difference. Stress terms are differenced in the conventional manner [82]. Equations 

(6.46) are used to advance the solution from time n to time n+1 and this process is 

continued until the desired integration time has been reached. 

The magnitude of the time step in Eqs. (6.46) is chosen based on the physical time 

scales present at any given time in the solution. The fluid dynamic time step, Atf, can 

be shown to be limited by the CFL condition [83] 

At < / Kr + ufj, | + 

/ - \ A{ Ari 

+ W 1 ♦ (m * if) 

where a is the local speed of sound. The chemical relaxation time for species i is given 
by [84] 


t c = 


Eli 

Wi 


(6.48) 
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Changes in this relaxation time are then given by 

A (pfi) 


A t r = 


wi 


(6.49) 


since w{ remains nearly constant over a time step. For accuracy, it is required that the 
chemical time step be chosen such that no change in mass fraction greater than 0.01 
occurs over that time step. The computational time step At is then chosen to be the 
minimum of all the grid points in terms of the both fluid and chemical time steps, i.e.. 


At = min (A fy, At c ) 


(6.50) 


6.2.3 Boundary and Initial Conditions 

The governing equations, Eqs. (6.1), require boundary conditions along all four 
boundaries. For the problems to be considered, the inflow boundary is supersonic, so the 
velocities, static temperature, pressure, and species mass fractions are specified and fixed 
there. The outflow boundary is also supersonic, and the values of the velocities, static 
temperature, pressure, and species mass fractions are determined by extrapolation from 
upstream values. Only the upper half of the flow domain is computed due to the assumed 
symmetry of the flow. The upper boundary is treated as a solid wall. This implies a non- 
slip boundary condition. The wall temperature is given and wall species mass fractions 
and pressures are extrapolated from interior grid points, by assuming a non-catalytic wall 
as well as the boundary layer assumptions on the pressure gradient. Symmetry boundary 
conditions are imposed at the lower boundary, that is, at the centerline. 

Equations (6.1) also require a set of initial conditions. The equations are initialized 
by setting values of the velocities, static temperature, pressure, and species mass fractions 
throughout the domain to the values chosen initially for boundary conditions at the inflow 
boundary. Having specified all required initial and boundary data, the equations are 
marched in time until steady state solutions are achieved. 
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6.2.4 Artificial Viscosity 

With the numerical algorithm of Sec. 6.2.2, high frequency nonlinear instabilities can 
appear as the solution develops. For example, flow oscillations can result from the odd- 
even decoupling inherent in the use of second-order central differencing for the inviscid 
terms. In addition, physical phenomena such as shock waves can cause instabilities when 
they are captured by the finite difference algorithm. Artificial viscosity, or smoothing, is 
normally added to the solution algorithm to suppress these high frequency instabilities. In 
this study, the artificial viscosity F av is added to the vector F in Eq. (6.41) as following 

Fav = (A [Crf\P + C 2 SjT + C^lf] (Uij - Ui-u) (6.51) 


where 


A = 


|{*ti + ( y v\ + ayJil + Q | Vxu + r] y v\ + + 

Af + At/ 


(6.52) 




P = 


\Pj+i,j 2Pjj + 

Pi+l,j 4* 2P,-J + Pi-1, j 


(6.53) 


* Ti+U - 2 Tij + Ti-ij 


(6.54) 



l/t+l J 2 fij + fi-l,j 
/*+!,> “ 2 fij + fi-l,j 


(6.55) 


The other artificial viscosity G av follows similar formulas as F av . Equation (6.51) was 
suggested by Pulliam [85]. In its original form, only the term was used; Singh et 
al. [86] found that for some problems especially those with chemical reaction this is not 
sufficient and suggested inclusion of the term S^T and In the term 6|/, f can be the 
mass fraction for one species or for several different species. The coefficient Ci, C2 and 
C3 must be selected by numerical experiment For the cases investigated in this study, 
all coefficients were fixed as a constant value of one half. 
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6.2.5 Solution Procedures 

With consideration of radiative heat transfer, solution procedures employed in this 
study are summarized as following: 

(a) First, Eqs. (6.1) were solved without consideration of radiation in terms of the 
modified MacCormack schemes; 

(b) The steady solutions for temperature, pressure and species mass fractions were 
then used in the Monte Carlo simulation. The computed radiative source term — V.(? r 
from the MCM was based on a different grid from that used for Eqs. (6.1). Linear 
interpolation and extrapolation were employed for the transformation of — V.g r between 
the two grids; 

(c) The transformed — V.g r was substituted into Eqs. (6.1), and Eqs. (6.1) were 
solved again. If the differences between two consecutive steady solutions were smaller 
than a designated tolerance, the computation was terminated. Otherwise, steps (b) and 
(c) were repeated until solutions converge. 

It is noted that there are two levels of numerical procedures employed here which 
result in two different iterative procedures. One is the numerical procedure for solving 
Eqs. (6.1) and their solutions were iterated with time. The other is the numerical 
procedure for evaluating the radiative source term using the MCM, which results in the 
iteration of steady state solutions. 

63 Results and Discussion 

Based on the theoretical and numerical analyses described earlier, a computer code 
has been developed to simulate two-dimensional supersonic chemically reacting and 
radiating nozzle flows on a Cray X-MP machine. The specific goal in this study was to 
investigate the effects of radiation on the flowfield and heat flux on the nozzle wall. By 
referring to [2], several problems have been considered. They contain four parameters: 
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equivalence ratio of hydrogen and air, inlet flow temperature, wall temperature and nozzle 
size. Numerical solutions have been obtained for a variety of combinations of these 
parameters. In each problem, flow is introduced at the nozzle as a uniform velocity of 
1230 m/s and a pressure of 1 atm. The grid size for solving the governing equations was 
71x41 (upper half of the nozzle). Further refinement of the grid produced only small 
changes in the results. For a given radiative source distribution, the residuals of Eqs. (1) 
were reduced by eight orders of magnitude in 3,000 iterations for a typical case and the 
steady state solutions were considered to have been obtained. The corresponding CPU 
time is about six minutes. 

To check the accuracy of the computational scheme, a preliminary calculation has 
been carried out for a chemically reacting nozzle flow without consideration of radiation 
and the present solution is compared with that in Ref. 87. Figure 6.4 shows the physical 
model for this calculation. It is noted that the nozzle walls are adiabatic walls in this case. 
Figures 6.5 and 6.6 demonstrate the frozen and reacting temperature distributions along 
the centerline. The present solution is found to agree with with the available solution [87]. 

For the temperature ranges considered, the important radiating species are OH and 
H 2 O. But OH is a much less radiation participating species compared to H 2 O. In addition, 
the concentration of OH is several times less than that of H 2 0 for the problems considered. 
So, the contribution of radiation from OH has been neglected in this study. For H 2 O, 
there are five important absorption bands. All these bands have been taken into account 
and they consist of 295 narrow bands in the spectral range from 150 cm -1 to 7500 cm -1 
[33]. In addition, for all the problems considered, the nozzle wall is assumed to be gray 
and the wall emissivity is taken to be 0.8. The inlet and outlet surfaces of the nozzle 
flow are treated as pseudoblack walls with the same temperatures as the local gases. 

To assure that the statistical results make sense in the Monte Carlo simulation, two 
requirements, the accuracy of the statistical results and the independence of the results on 
the grid, must be met. In this study, the designated statistical accuracy of the results is 
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Fig. 6.4 Physical model for validation calculation. . 
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Fig. 6.5 Comparison of frozen temperatures along the centerline. 
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Fig. 6.6 Comparison of reacting temperatures along the centerline. 
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defined in such a way that when the relative statistical errors of results are less than ±5%, 
the probability of the results lying within these limits is greater than 95%. Independence 
of the results on a grid is considered to have been achieved when the volume element 
number in the x direction is 20 and the volume element number in the y direction is 20 as 
shown in Fig. 6.3. For this grid, the total number of energy bundles had to be 5,000,000 
and the required CPU time was about one hour in order to meet the designated statistical 
accuracy in results for a typical problem. To test the independence of the Monte Carlo 
results on the grid, the same problem was investigated with a finer gird in which the 
volume element number in the x direction was increased to 30 and the volume element 
number in the y direction was doubled. To obtain the same accurate results, the total 
number of energy bundles had to increase to 15,000,000 and the corresponding CPU 
time increased to three hours. Comparing the solutions for the two different grids, it 
was found that the difference for the net radiative wall flux was never more than 2%, 
and the difference for the radiative source term was a little higher but less than 10%. In 
fact, the net radiative wall flux is the quantity we are most interested in, and its accuracy 
seems more important to us. 

The grid considered for Monte Carlo computations in this study is coarser than that 
for numerical solutions of the energy equation. The intermediate values of the radiative 
source term within the grid for solutions of Eqs. (6.1) are obtained by interpolation 
and extrapolation. This should not introduce significant errors as the radiative source 
term is a slowly varying function compared to the temperature and its derivatives [6]. 
The major CPU time consumed is in the Monte Carlo simulation. Fortunately, Monte 
Carlo subroutines only need to be called one or two times to obtain converged steady state 
solutions. The reason for this will be explained later. It is believed that the computational 
time for Monte Carlo simulations can be reduced considerably if the code is vecterized 
and parallelized. 
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The radiative effects on the flowfield were investigated first It is common knowledge 
that convective heat transfer is very strong for a supersonic flow. Hence the effects 
of radiation may not be very important. To determine these effects quantitatively, a 
typical problem was selected in which the equivalence ratio of hydrogen and air, wall 
temperature, inlet flow temperature and the nozzle length are taken to be 0=1.0, T w =1900 
K, Tj=1900 K and L=2.0 m. The inlet species mass fractions are fjj 2 = 0.0283, fo 2 — 
0.2264, S H2 o = 0.0, }oh = 0.0, fo = 0.0, f H = 0.0, f N2 = 0.74529. Figures 
6.7-6.10 show the temperature, pressure, density distributions and velocity vector plots, 
respectively. Figures 6.11-6.14 show the mass fraction distributions for the species H 2 O, 
OH, H 2 and O 2 , respectively. Knowing this information is essential in analyzing the 
effect of radiative heat transfer. As the premixed mixture of hydrogen and air enters the 
nozzle, an exothermic chemical reaction takes place immediately, and the temperature, 
pressure, and density increase abruptly, reaching their peaks in a region closer to the 
inlet location (Figs. 6.7-6.9) while the velocity decreases slightly (Fig. 6.10). During 
this rapid change in temperature, pressure, and density, the two major products H 2 O and 
OH experience a big jump in mass fraction (Figs. 6.11 and 6.12) while the two major 
reactants H 2 and O 2 experience a big drop in mass fraction (Figs. 6.13 and 6.14). As 
the flow continues to move downstream, supersonic expansion plays a major role, and 
the temperature, pressure, and density are decreased while velocity is increased. At the 
same time, the chemical reaction proceeds but it becomes very weak. This is why there 
is little change in mass fractions for the species H 2 0, OH, H 2 and O 2 in the downstream 
region. Computation has been conducted for other cases also. Similar trends in results 
for temperature, pressure, density, velocity, and mass fractions for all species were also 
observed. 

Figure 6.15 shows the radiative source distributions at three different locations for 
the case considered in Figs. 6.7-6.14. At the location x/L=0.1, temperature and pressure 
are very high and there is more radiant energy emitted than absorbed. Consequently, the 
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radiative source distribution is higher than at locations x/L=0.5 and 0.9. The trend in 
results for — V.<? r at the location x/L=0.1 is seen to be different from the results of other 
locations due to a decrease in temperature as the distance from the center line increases. 
The convective heat transfer distributions for the same locations as in Fig. 6.15 were 
also calculated but they are not plotted in Fig. 6.15. This is because of large differences 
between the convective and radiative results; and also due to opposite signs for convective 
results at different locations. In most regions, the absolute value of the convective heat 
transfer is two or three orders of magnitude larger than the radiative source term. This 
situation does not change as long as the speed of the flow is very high. Consequently, the 
effects of radiation on the flowfield are very weak for supersonic flows. This confirms 
our expectation and also answers the question that the Monte Carlo subroutine only needs 
to be called one or two times to obtain converged steady state solutions. As a matter of 
fact, a case without radiation was considered and the differences in temperature, pressure, 
density, and species mass fractions between the two cases were found to be less than ±1%. 

The radiative effects on the heat transfer on the nozzle walls are investigated next. 
Unlike the radiative effects on the flowfield, the effects of radiation on the nozzle wall 
flux are significant when compared with those from conduction. The following results 
will demonstrate the relative importance of radiative and conductive wall fluxes and how 
they change with equivalence ratio, wall temperature, inlet flow temperature, and nozzle 
size. Here, the conductive wall flux is defined as 

)_ «« 

where n represents normal direction to the wall. 

The effects of the equivalence ratio <^> on q rw and q cw are illustrated in Fig. 6.16. For 
a specific <f> value, q cw is seen to increase first, reach a peak and then decrease. This is 
compatible with the trends in temperature variation as seen in Fig. 6.7. Unlike q cw , q rw 
is seen to increase with distance along the nozzle. This behavior is justifiable. In this 
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study, the inlet and outlet of the flow are treated as the pseudoblack walls. The outlet flow 
temperatures are larger than the inlet flow temperatures and the outlet area is also bigger 
than the inlet area. In addition, as the flow goes downstream, the cross-sectional area of 
the flow increases. Consequently, the optical length increases. These two reasons result 
in higher value of q rw as the distance from the inlet location increases. Comparing the 
values of q TW and q cw for each case, it is clear that radiation is predominant. Even in the 
inlet region, q rw is more than two times higher than q cw . The results for three different 
equivalence ratios reveal different behavior for combustion with lean and rich mixtures. 
As (j> increases from 0.6 to 1.0, the flow temperature and H 2 O mass fraction increase by 
about 10% and 50% respectively, and pressure decreases by about 5%. The effects of 
these changes result in a sizable increase in the values of q rw and q cw . However, as <j) 
increases from 1.0 to 1.4, the flow pressure decreases by about 5% and the H 2 O mass 
fraction increases by about 15%, but the temperature shows little change. This results in 
only a slight change in the values of q rw and q cw . 

Figure 6.17 shows the effects of the nozzle wall temperature on q rw and q cw . The 
change of the nozzle wall temperature is found to have little influence on the combustion, 
and the flow temperature, pressure and H 2 O mass fraction remain almost the same in 
most regions as T w varies from 1500 K to 2100 K. As a result, the magnitude of the 
radiant energy absorbed on the wall is very close for the three cases with different nozzle 
wall temperatures. The value of q rw is equal to the absorbed radiant energy minus the 
emitted radiant energy. So q rw is reduced with higher wall temperature, as seen in Fig. 

6.17. As far as q cw is concerned, except in the entrance region , q cw is seen to exhibit 
minor changes among the cases with different wall temperatures. 

The effects of the inlet flow temperature on q rw and q cw are demonstrated in Fig. 

6.18. Inspection of the distribution of the q rw value among the three cases reveals a 
very interesting feature of q TW . The values of q rw along the wall are not monotonically 
increased or decreased with Ti. The combined effects of temperature, pressure and H 2 O 
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mass fraction in the flow on radiation are responsible for this behavior. It is well known 
that increase of temperature, pressure and concentration of participating media enhances 
radiation. As Tj varies from 1500 K to 1800 K and then from 1800 K to 2100 K, the flow 
temperature increases by about 5% while the pressure and H 2 O mass fraction decrease 
by about 10% and 15% respectively at each stage. An increase in temperature tries to 
reinforce the radiation while a decrease of pressure and H 2 O mass fraction tries to reduce 
the radiation. So there exist two driving forces which compete with each other to affect 
the radiation. As a consequence of the competition, the lowest curve for q rw is seen for 
the case with T t = 1800 K and the highest values are observed for the case with Ti= 1500 
K. Unlike q rw , the values for q cw are found to increase monotonically with TV This is 
because the convective wall flux is only dependent on temperature. 

Finally, the effects of the nozzle size on q rw and q cw are illustrated in Fig. 6.19. 
By changing the nozzle length, geometrically similar nozzles with different sizes can be 
obtained. As the nozzle length is reduced from 2.0 m to 1.0 m and then from 1.0 m to 
0.5 m, the flow temperature and H 2 O mass fraction are decreased by about 5% while the 
pressure is increased by about 2% at each stage. The effect of increased pressure on the 
radiation is overshadowed by a decrease in the nozzle size, temperature and H 2 O mass 
fraction. Hence, lower values of q rw are seen in the figure as the nozzle length is reduced. 
For the smaller nozzle size, the flow temperature may be lower, but the normal derivative 
of temperature is actually higher. Therefore, contrary to q rw , the value q cw is observed 
to increase with a decrease in the nozzle size. The opposite trend between the values of 
q rw and q cw brings a question about the role of radiation in heat transfer on the nozzle 
wall. With a decrease of nozzle size, the differences between the values of q rw and Qcw 
are reduced and the dominance of radiation is diminished. In fact, at L=0.5, the value 
of q cw is larger than the value of q rw in some parts of the nozzle wall. It is expected 
that radiation will become less important and conduction will replace radiation as the 
dominant mode of heat transfer on the nozzle wall if the nozzle size is reduced further. 
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Chapter 7 

CONCLUDING REMARKS 

The MCM has been applied to investigate radiative heat transfer in a nongray 
participating medium in an exact manner. The nongray model employed is based on 
a random statistical narrow band model. When a narrow band model is employed in the 
MCM, the spectral correlation only occurs between the transmittances of two different 
segments of the same path in the statistical relationship for determining the absorption 
location of an energy bundle. For the case with reflecting walls, Monte Carlo treatment 
with a narrow band model is similar to that with a gray model, and the spectral correlation 
between the reflected component of the wall radiosity and the transmittance occurring in 
other methods does not exist. Consideration of different problems reveals that the Monte 
Carlo solutions are in good agreement with available results of other methods but the 
MCM is much simpler to implement than other methods. 

The validity of the Monte Carlo correlated formulations is further established by 
considering the steady-state energy transfer in laminar, incompressible, constant proper- 
ties, fully developed flow of absorbing-emitting gases between two parallel plates. The 
nongray Monte Carlo solutions were found to be in good agreement with the available 
approximate solutions. The gray Monte Carlo solutions were also obtained for the same 
problem and they also essentially match the available analytical solutions. 

The exact correlated and non-correlated Monte Carlo formulations are very com- 
plicated for multi-dimensional systems. Solutions of these formulations are extremely 
difficult, if not impossible. However, by introducing the assumption of an infinitesi- 
mal volume element, the approximate correlated and non-correlated formulations were 
obtained which were tractable compared to the exact formulations. Consideration of dif- 
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ferent problems and comparison of different solutions reveals that the approximate and 
exact correlated solutions agree very well, and so do the approximate and exact non- 
conelated solutions. However, the two non-correlated solutions lack physical meanings 
because they usually differ from the correlated solutions significantly. An accurate pre- 
diction of radiative heat transfer in any nongray and multi-dimensional system is possible 
by using the approximate correlated formulations. 

By investigating the radiative interactions for chemically reacting supersonic flows of 
premixed hydrogen and air in an expanding nozzle, the correlated Monte Carlo method 
developed earlier has been found to be a very convenient and reliable tool to analyze 
radiative heat transfer in multi-dimensional nongray systems. For chemically reacting 
supersonic flows, the effects of radiation on the flowfield can be neglected but the radiative 
effects on the heat transfer on the nozzle wall are significant The extensive parametric 
studies on the radiative and conductive wall fluxes have demonstrated that the magnitude 
of the radiative and conductive wall fluxes are very sensitive to the equivalence ratio when 
the equivalence ratio is less than 1.0 but they are less sensitive when the equivalence ratio 
is higher than 1.0. The change in the wall temperature has little effect on the combustion. 
Thus, the radiative wall flux is decreased with increases in wall temperature. But the 
conductive wall flux seems insensitive to changes in wall temperature. The radiative 
wall flux does not change monotonically with inlet flow temperature. Lower inlet flow 
temperature can yield higher radiative wall flux. The conductive wall flux, however, 
increases with an increase in the inlet flow temperature. The radiative wall flux decreases 
but the conductive wall flux increases with a reduction in nozzle size. For larger nozzles, 
the radiative wall flux is dominant over the conductive wall flux. However, the situation 
can be reversed when the nozzle size is reduced. 
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APPENDIX A: 

PROGRAM LISTING FOR MONTE CARLO SIMULATION 

This code is developed to calculate radiative source distribution and net radiative 
wall flux in a one-dimensional problem. A lot of subroutines in the code are from the 
IMSL library. The spectral correlation has been taken into account. Input files consist 
of parameter statement file “param.inc”, common statement file “common.inc”, and three 
narrow band information files “y*dat”, “fibig.dat”, and “f2big.dat”. Temperature, pressure, 
and concentration distributions should be also given before calculation, 
program moncar 
include ’paramm.inc’ 
include ’commonm.inc’ 
parameter (mx=22,mz=22) 
external gamfun,bs2vl,funtao 
real rwksp(20000),tarray(2),len(3) 
dimension rf(mx),t2(mx),xl(mx) 
dimension sg(mx),sq(mx),nn(mx),em(3) 
dimension xp(mx),zp(mz),tm(mx,mz) 
common/cgas/p,xh2o,xn2,xo2,xco2,dlx 
common/worksp/ rwksp 
common/ct/tp 

data xh2o,xn2,xo2,xco2/l .0,0.0, 0.0,0. 0/ 
data em/0.90,0.8,0.1/ 


data len/1.0, 60., 100.0/ 



133 


c open(unit=4,file="tinbig.dat") 
open(unit=5,file="tout.dat") 
call iwkin(20000) 
call coefbs 

c read(4,*) (t2(i),i=l,mx) 
dlx= 1 .0/float(mx-2) 
dlz=l .0/float(mz-2) 
xp(l)=0.0 
xp(2)=0.5*dlx 
xp(mx)=1.0 
zp(l)=0.0 
zp(2)=0.5*dlz 
zp(mz)=1.0 
do 1 i=3,mx-l 

1 xp(i)=xp(i-l)+dlx 
do 2 j=3,mz-l 

2 zp(j)=zp(j-l)+dlz 
do 3 i=l,mx 

do 3 j=l,mz 

tm(i,j)=500.0+500.0*( 1 .0-(2.0*zp(j)- 1 .0)**2)*xp(i) 

3 continue 

do 303 j=l,mz 
tm(mxj)=300.0 
303 continue 


do 999 i9=3,mx,2 
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do 304 j=l,mz 
t2(j)=tm(i9,j) 

304 continue 

time 1 =etime(tarray) 
pi=3. 14 15926 
n=100000 

p=1.0 

do 910 il=2,2 
el=em(il) 
e2=em(il) 
do 900 i2=2,2 
alx=len(i2) 

dlx=alx/(real(mx)-2.0) 

nran= 15249649 

call mset (nran) 

sum =0.0 

do 20 i=l,mx 

rf(i)=0.0 

sg(i)=0. 

20 sq(i)=0.0 

xl(2)=0.5/(real(mx)-2.0) 
do 5 i=2,mx-l 

xl(i)=xl(2)+float(i-2)/(real(mx)-2.0) 

t=t2(i) 

call baneng(t,l) 



sq(i)=qv 

sum=sum+sq(i) 

5 continue 
t=t2(l) 
call surwc(t) 
sq(l)=sqw*el 
sq(mx)=sqw*e2 
sum=sum+sq(l)+sq(mx) 
sqm=sum/float(n) 
do 13 i=l,rax 
nn(i)=ifix(sq(i)/sqm+0.5) 
13 continue 
isl=l 
ntt=0 
nt=0 

23 go to (30, 3 1,36, 120), is 1 

30 it=l 

go to 39 

31 it=mx 

go to 39 
36 it=2 
38 t=t2(it) 

gamma=gamfun(t) 
call baneng(t,2) 


39 is2=is 1 



ar=munf() 

go to (40,4 l,45),isl 

40 ran=munf() 

ran=ws 1 +ran*(ws2-ws 1 ) 
waveno=bsval(ran,kxord,xks 1 ,nxd 1 ,bss 1 ) 
go to 471 

41 ran=munf() 

ran=ws 1 +ran*(ws2- ws 1 ) 
waveno=bsval(ran,kxord,xks 1 ,nxd l,bss 1 ) 
go to 471 
45 ran=munf() 

waveno=bsval(ran,kxord,xkv 1 ,nxd 1 ,bsv 1 ) 

47 park=bs2vl(waveno,t,kxord,kyord,xkl ,ykl ,nxd 1, 
+nydl,bscoll) 

pardlt=bs2vl(waveno,t,kxord,kyord,xk 1 ,yk 1 ,nxd 1 , 
+nydl,bscol2) 

471 nt=nt+l 
i=it 

go to (48,48,49),isl 

48 u=0.0 
sumk=0.0 
sumb=0.0 
go to 50 

49 ranl=munf() 


if(ranl.gt.0.5) go to 491 



rx=bsval(ran 1 ,kxord,xkb 1 ,nb,bsb 1) 
go to 492 

491 ranl=1.0-ranl 

rx=-bsval(ran 1 ,kxord,xkb 1 ,nb,bsb 1 ) 

492 rxl=dlx/abs(rx) 
u=xh2o*p*rxl 
beta=2.0*gamma*pardlt 
ar=(l .0-funtao(u,park,beta))*ar 
ul=0.0 

sumk 1=0.0 
sumb 1=0.0 
us=u 

sumks=park*u 
sumbs=sumks*beta 
if(rx.lt.O.) go to 496 
494 i=i+l 

if(i.gt.(mx-l)) go to 80 
ta=t2(i) 

gama=gamfun(ta) 

park=bs2vl(waveno,ta,kxord,kyord,xk 1 ,yk 1 ,nxd 1 , 
+nydl,bscoll) 

pardlt=bs2vl(waveno,ta,kxordJcyord,xkl ,ykl ,nxd 1 , 

+nydl,bscol2) 

deltu=p*rxl*xh2o 


ul=ul-i-deltu 



sumkl=sumkl+park*deltu 

sumbl=sumbl+park*deltu*2.0*gama*pardlt 

u2=ul+us 

sumk2=sumkl +sumks 
sumb2=sumb 1 +sumbs 
efkl=sumkl/ul 
efb 1 =sumb 1/u 1/efk 1 
efk2=sumk2/u2 
efb2=sumb2/u2/efk2 

dtao=funtao(u 1 ,efkl ,efbl)-funtao(u2,efk2,efb2) 
if(dtao.lt.ar) go to 85 
go to 494 
496 i=i-l 

if(i.lt.2) go to 79 
ta=t2(i) 

gama=gamfun(ta) 

park=bs2vl(waveno T ta,kxord,kyord,xk 1 ,yk 1 ,nxd 1 , 
+nydl,bscoll) 

pardlt=bs2vl(waveno,ta,kxord,ky ord,xk 1 ,yk 1 ,nxd 1 , 

+nydl,bscol2) 

deltu=p*rxl *xh2o 

ul=ul+deltu 

sumk 1 =sumk 1 +park*deltu 

sumb 1 =sumb l+park*deltu*2.0*gama*pardlt 


u2=ul+us 



sumk2=sumk 1 +sumks 


sumb2=sumb 1 +sumbs 
efkl=sumkl/ul 
efb 1 =sumb 1/u 1/efkl 
efk2=sumk2/u2 
efb2=sumb2/u2/efk2 

dtao=funtao(u 1 ,efkl ,efb l)-funtao(u2,efk2,efb2) 
if(dtao.lt.ar) go to 85 
go to 496 
50 ranl=munf() 
go to (54,55), is2 
54 rx=sqrt(l.-ranl) 

541 i=i+l 

if(i.gt.(mx-l)) go to 80 
ta=t2(i) 

gama=gamfun(ta) 

park=bs2vl(waveno,ta,kxord,kyord,xk 1 ,yk 1 ,nxd 1 , 
+nydl,bscoll) 

pardlt=bs2vl(waveno,ta,kxord,ky ord,xk 1 ,yk 1 ,nxd 1 , 
+nydl,bscol2) 
deltu=p*xh2o *dlx/rx 
go to (543,543,544),isl 
543 u=u+deltu 

sumk=sumk+park*deltu 

sumb=sumb+park*deltu s|e 2.0*gama*pardlt 



efk=sumk/u 


efb=sumb/u/efk 
tao=funtao(u,efk,efb) 
if(tao.lt.ar) go to 85 
go to 541 
544 ul=ul+deltu 

sumk 1 =sumk 1 +park*deltu 

sumbl=surnbl+park*deltu*2.0*gama*pardlt 

u2=ul+us 

sumk2=sumk 1 +sumks 

sumb2=sumbl+sumbs 

efkl=sumkl/ul 

efb 1 =sumb 1/u 1/efk 1 

efk2=sumk2/u2 

efb2=sumb2/u2/efk2 

dtao=funtao(ul,efkl,efbl)-funtao(u2,efk2,efb2) 
if(dtao.lt.ar) go to 85 
go to 541 
55 rx=-sqrt(l.-ranl) 

551 i=i-l 

if(i.lt.2) go to 79 
ta=t2(i) 

gama=gamfun(ta) 

park=bs2vl(waveno,ta,kxord,kyord,xk 1 ,yk 1 ,nxd 1 , 
+nydl,bscoll) 
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pardlt=bs2vl(waveno,ta,kxord,kyord,xk 1 ,yk 1 ,nxd 1 , 

+nydl,bscol2) 

deltu=-p*xh2o*dlx/rx 

go to (553,553,554), is 1 

553 u=u+deltu 
sumk=sumk+park*deltu 
sumb=sumb+park*deltu*2.0*gama*pardlt 
efk=sumk/u 

efb=sumb/u/efk 
tao=funtao(u,efk,efb) 
if(tao.lt.ar) go to 85 
go to 551 

554 ul=ul+deltu 

sumk 1 =sumk 1 +park*deltu 

sumb 1 =sumb 1 +park*deltu*2.0*gama*pardlt 

u2=ul+us 

sumk2=sumk 1 +sumks 
sumb2=sumb 1 +sumbs 
efkl=sumkl/ul 
efb 1 =sumb 1/u 1/efk 1 
efk2=sumk2/u2 
efb2=sumb2/u2/efk2 

dtao=funtao(u 1 ,efk 1 ,efb l)-funtao(u2,efk2,efb2) 
if(dtao.lt.ar) go to 85 
go to 551 
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79 is2=l 

go to 86 

80 is2=2 

go to 86 

85 is2=3 

86 go to (90,91,1 10) ,is2 

90 i=l 

go to 92 

91 i=mx 

92 ran=munf() 

if((is2.eq.l.and.ran.gt.el).or.(is2.eq.2.and.ran.gt.e2)) 
+go to 50 
sg(i)=sg(i)+sqm 
59 if(nt.lt.abs(nn(it))) go to 39 
ntt=ntt+nt 
nt=0 

102 if(isl.ne.3) go to 104 
it=it+l 

if(it.le.(mx-l)) go to 38 
104 isl=isl+l 
go to 23 

110 sg(i)=sg(i)+sqm 
go to 59 

120 call sub200(sg,mx) 


qw 1 =-(sq( 1 )-sg( 1 ))*0.00 1 
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qw2=-(sq(mx)-sg(mx))*0.001 
write(5,25) el,e2,alx,n 

25 format(lx/el=\f63,2x/e2=\f6.3,2x,’x=\f9.3,2x,’n=’,i9/) 
write(5,26) qwl,qw2 

26 format(lx/qwl==\fl5.6,4x,’qw2==’,fl5.6/) 
do 130 i=2,mx-l 
rf(i)=-sg(i)+sq(i) 
rf(i)=-rf(i)/dlx*0. 1 

write(5,27) xl(i),rf(i) 

27 format(lx,fl5.6,4x,fl5.6) 

130 continue 

time2=etime(tarray) 
time=time2-timel 
write(5,390) time 

390 format(lx,’cpu time spent =\f9.3 ////) 

900 continue 
910 continue 
999 continue 
stop 
end 
c 
c 

subroutine coefbs 


external bs2in,bsnak 
include ’paramm.inc’ 



include ’commonm.inc’ 


dimension fd 1 1 (nxd 1 ,nyd 1 ),fd 1 2(nxd 1 ,nyd 1 ),yd 1 (nyd 1 ) 

open(unit=7,file="y.dat") 

open(unit=8,file="f lbig.dat") 

open(unit=9,file="f2big.dat") 

read(7,*) (yd 1 (i),i= 1 ,nyd 1 ) 

read(8 , * ) (i,xd 1 (i),(fd 1 1 (i,j ),j = 1 ,nyd 1 ) ,i= 1 ,nxd 1 ) 

read(9,*) (i,xd 1 (i),(fd 1 2(i,j ),j= 1 ,nyd 1 ),i= 1 ,nxd 1 ) 

call bsnak(nxdl,xdl,kxord,xkl) 

call bsnak(nydl,ydl,kyord,ykl) 

call bs2in(nxd 1 ,xd 1 ,nyd 1 ,yd 1 ,fd 1 1 ,ldf 1 dcxord, 

+kyord,xkl,ykl,bscol 1) 

call bs2in(nxd 1 ,xd 1 ,nyd 1 ,yd 1 ,fd 1 2,ldf 1 ,kxord, 

+ky ord,xk 1 ,yk 1 ,bsco 1 2) 
return 
end 
c 
c 

subroutine baneng(t,iflag) 

include ’paramm.inc’ 

include ’commonm.inc’ 

external planck,gamfun,bs2vl,emicoe,bsint 

dimension bre 1 (nxd 1 ),cpmu 1 (nxd 1 ,nb) 

dimension tgl(nxdl) 

common/cgas/p,xh2o,xn2,xo2,xco2,dlx 



common/cpar/gamma,park,pardlt 

common/ct/tp 

eps=1.0e-03 

gamma=gamfun(t) 

tp— t 

delomg=xd 1 (2)-xd 1(1) 
dmu=1.0/float(nb-l) 
do 5 i=l,nb 

mu(i)=l .0-dmu*float(i- 1) 

5 continue 
do 30 i=l,nxdl 
x=xdl(i) 
r=planck(x) 

park=bs2vl(x,t,kxord,kyord,xkl ,ykl ,nxd 1 , 
+nydl,bscoll) 

pardlt=bs2vl(x,t,kxord,kyord,xkl ,ykl ,nxd 1 , 
+nydl,bscol2) 
go to (25,10) iflag 
10 do 15 j=l,nb 
xa=mu(j) 

call qdags(emicoe,1.0,xa,eps,eps,rl,err) 
cpmu 1 (i,j )=2.0*r 1 *r*delom g 
15 continue 

25 call qdags(emicoe,1.0,0.0,eps,eps,rl,err) 


bre 1 (i)=4.0*rl *r*delomg 



30 continue 


qv=0.0 

do 40 i=l,nxdl 
qv=qv+brel(i) 

40 continue 

go to (50,55) iflag 
50 go to 99 
55 eps4=1.0e-7 
d=brel(l) 
tgl(l)=0.0 
do 60 i=2,nxdl 
d=d+brel(i) 
tgl(i)=d/qv 

if((tg 1 (i)-tg 1 (i-l)).le.eps4) tg 1 (i)=tg 1 (i- 1 )+eps4 
60 continue 

call bsnak(nxdl,tgl,kxord,xkvl) 

call bsint(nxd 1 ,tg 1 ,xd 1 ,kxord,xkv 1 ,bs v 1 ) 

do 70 j=l,nb 

cpmu(j)=0.0 

do 80 i=l,nxdl 

cpmu(j)=cpmu(j)+cpmu 1 (i,j) 

80 continue 

cpmu(j)=cpmu(j)/qv 
70 continue 


call bsnak(nb,cpmu,kxord,xkbl) 


call bsint(nb,cpmu,mu,kxord,xkb 1 ,bsb 1) 

99 return 
end 
c 
c 

function emicoe(x) 

common/cgas/p,xh2o,xn2,xo2,xco2,dlx 

comm on/cpar/gamma,park,pardlt 

eti=sqrt( 1 .0+xh2o*p*dlx*park/gamma/pardlt/x)- 1 .0 

tao=exp(-2.0*gamma*pardlt*eti) 

emicoe=-( 1 .0-tao)*x 

return 

end 

function gamfun(t) 

common/cgas/p,xh2o,xn2,xo2,xco2,dlx 

ts=296.0 

gamfun=0.066*p*(7.0*xh2o*ts/t+sqrt(ts/t)*( 1 .2*(xh2o+xn2)+ 
+0.8*xo2+ 1 .6*xco2)) 
return 
end 
c 
c 

subroutine surwc(t) 
include ’paramm.inc’ 


include ’commonm.inc’ 



external planck,bsint,qdags 
dimension tl(nxdl) 
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common/ct/tp 

eps=1.0e-04 

x 1=99999.0 

s=5.6696e-08 

power=s*t**4 

tp=t 

call qdags(planck,0.0,xdl(l),eps,eps,rl,err) 

call qdags(planck,0.0,xd 1 (nxd 1 ),eps,eps,r2,err) 

sqw=r2-rl 

wsl=rl /power 

ws2=r2/power 

eps4=1.0e-7 

tl(l)=wsl 

do 10 i=2,nxdl 

x=xdl(i) 

call qdags(planck,0.0,x,eps,eps,re,err) 
tl(i)=re/power 

if(t 1 (i).le.tl(i- 1 )) tl (i)=t 1 (i- 1 )+eps4 
10 continue 

call bsnak(nxdl,tl4cxord,xksl) 

call bsint(nxd 1 ,t 1 ,xd 1 ,kxord,xks 1 ,bss 1 ) 
return 


end 
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C 

c 

function planck(x) 
common/ct/tp 
cl=3.740e-08 
c2=1.4387 

planck=c 1 *x* *3/ (exp(c2*x/tp)- 1 .0) 

return 

end 

function funtao(u,park,beta) 
funtao=exp(-beta*(sqrt( 1 .0+2.0*u*park/beta)- 1 .0)) 
return 
end 
c 
c 

subroutine sub200(sg,mx) 
dimension sg(990) 
ms=(mx+l)/2 
do 210 i=l,ms 
is=mx-i+l 

sg(i)=(sg(i)+sg(is))/2. 

210 sg(is)=sg(i) 

240 return 


end 
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APPENDIX B: 

PROGRAM LISTING FOR RADIATIVE INTERACTIONS 
IN LAMINAR FLOWS USING MONTE CARLO SIMULATION 

This code is developed to investigate the radiative interactions in laminar flows 
between two parallel plates. The Monte Carlo simulation subroutine “moncar” used 
here is the same as given in Appendix A. 
program mcom 
include ’param.inc’ 
include ’common.inc’ 
parameter (mx=100,n=20) 
real rwksp(20000),tarray(2) 
external fcn,neqnf 
dimension q(n),tn(mx),te(mx),qgr(n) 
dimension ald(10),qtd(3),temper(4),pres(4) 
common /al/qly,aky,h,alx,tw,q,qgr 
common/worksp/ rwksp 
common/cgas/p,xh2o,xn2,xo2,xco2,dlx 
data xh2o,xn2,xo2,xco2/l .0,0.0,0.0,0.0/ 
data w/0.5/,eps/ 1 .0e-04/ 
data pres/1.0,1.0,5.0,10.0/ 
data temper/300.0,500.0,1000.0,2000.0/ 
data qtd/1.0e06,5.0e05,1.0e07/ 
data ald/0.0 1 ,0.05,0. 1 ,0.5,1 .0,5.0, 10.0,20.0,50.0, 100.0/ 



open(unit=5,file="ngput.dat") 

call iwkin (20000) 

call coefbs 

n 1=50000 

errrel=1.0e-03 

itmax=50 

e=1.0 

time 1 =etime(tarray) 
do 599 ii=l,l 
p=pres(ii) 
do 690 ij=2,2 
tw=temper(ij) 

fb = (4 1 86.8/360* 130.)*((tw/273.0)** 1.48) 

do 691 it=2,2 

qt=qtd(it) 

do 692 kk=6,6 

alx=ald(kk) 

h=alx/(real(2*n)- 1 .0) 

dlx=h 

qly= 1 2.0*qt*h/( 1 000.0*alx) 

aky=fb/h/l 000.0 

qrt=qt*alx/fb 

iter=0 

sumt=0.0 


sumt 1=0.0 
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write(5,l) nl,n,alx,tw,p,qt 
1 format(//lx, , nl=\ilO,lx,’n=’,i3,2x,’alx=’,f9.3, 
*2x,lx/tw=\f9.3,2x,*p=\f6.2,2x,’qt=\lpe9.3) 
do 3 i=l,n 
tn(i)=tw 

3 continue 
tn(l)=tw 

99 iter=iter+l 

if(iter.gt.21) go to 98 

tiraea=etime(tarray) 

call moncar(nl,n,alx,e,tn,q) 

timeb=etime(tarray) 

delt=timeb-timea 

sumt=sumt+delt 

do 4 i=l,n 

te(i)=tn(i) 

call baneng(te(i),l) 
qgr(i)=qv 

4 continue 

timea 1 =etime(tarray ) 

call neqnf(fcn,errrel,n,itmax,te,tn,fnorm) 

timeb 1 =etime(tarray) 

delt 1 =timeb 1 -timea 1 

sumt 1 =sumt 1 +deltl 


do 6 i=l,n 



tn(i)=( 1 .0-w)*te(i)+w*tn(i) 

6 continue 
do 7 i=l,n 

if(abs((tn(i)-te(i))/tn(i)).gt.eps) go to 99 

7 continue 

98 write(5,ll) iter 

11 format(2x,’iter=\i5/) 
write(5,12) (i,tn(i),i=l,n) 

12 format(2x,i6,2x,lpel2.5) 
mi=2*n 

do 13 i=n+l,mi 
tn(i)=tn(mi-i+l) 

13 continue 

do 15 i=l,mi 

x=real(i-l)*h/alx 

tn(i)=(tn(i)-tw)/qrt*(x-x*x)*6. 

15 continue 
sum=0.0 
do 17 i=2,mi 
su=(tn(i)+tn(i- l))*h/2./alx 
sum=sum+su 

17 continue 
write(5,18) sum 

18 format^ bulk temperature for monte carlo solution=’ , 1 pe 1 1 A! Ill) 


time2=etime(tarray) 
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time=time2-timel 

timel=time2 

write(5,39) time,sumt,sumtl 
39 format(lx,’cpu time spent =\f9.3/lx, 

*’cpu time spent for monte carlo simulation=\f9.3/lx, 
*’cpu time spent for solving set of equations=\f9.3///) 
692 continue 
691 continue 
690 continue 
599 continue 
stop 
end 
c 
c 

subroutine fcn(tn,f,n) 

parameter (m=20) 

dimension tn(n),f(n),q(m),qgr(m) 

common /al/qly,aky,h,alx,tw,q,qgr 

tn(l)=tw 

f(l)=0.0 

do 2 i=2,n-l 

x=real(i-l)*h/alx 

f(i)=qgr(i)-(tn(i-l)-2.*tn(i)+tn(i+l))*aky 

*+qly*(x-x**2)-q(i) 


2 continue 
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x=real(n-l)*h/alx 

f(n)=qgr(n)-(tn(n- 1 )-2. * tn(n)+tn(n)) *aky 
*+qly*(x-x**2)-q(n) 
return 


end 















